source: ThirdParty/levmar/matlab/levmar.c@ 1ca493a

Candidate_v1.7.0 stable
Last change on this file since 1ca493a was 8ce1a9, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '5443b10a06f0c125d0ae0500abb09901fda9666b' as 'ThirdParty/levmar'
  • Property mode set to 100644
File size: 23.3 KB
Line 
1/* ////////////////////////////////////////////////////////////////////////////////
2//
3// Matlab MEX file for the Levenberg - Marquardt minimization algorithm
4// Copyright (C) 2007 Manolis Lourakis (lourakis at ics forth gr)
5// Institute of Computer Science, Foundation for Research & Technology - Hellas
6// Heraklion, Crete, Greece.
7//
8// This program is free software; you can redistribute it and/or modify
9// it under the terms of the GNU General Public License as published by
10// the Free Software Foundation; either version 2 of the License, or
11// (at your option) any later version.
12//
13// This program is distributed in the hope that it will be useful,
14// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16// GNU General Public License for more details.
17//
18//////////////////////////////////////////////////////////////////////////////// */
19
20#include <stdio.h>
21#include <stdlib.h>
22#include <stdarg.h>
23#include <math.h>
24#include <string.h>
25#include <ctype.h>
26
27#include <levmar.h>
28
29#include <mex.h>
30
31/**
32#define DEBUG
33**/
34
35#ifndef HAVE_LAPACK
36#ifdef _MSC_VER
37#pragma message("LAPACK not available, certain functionalities cannot be compiled!")
38#else
39#warning LAPACK not available, certain functionalities cannot be compiled
40#endif /* _MSC_VER */
41#endif /* HAVE_LAPACK */
42
43#define __MAX__(A, B) ((A)>=(B)? (A) : (B))
44
45#define MIN_UNCONSTRAINED 0
46#define MIN_CONSTRAINED_BC 1
47#define MIN_CONSTRAINED_LEC 2
48#define MIN_CONSTRAINED_BLEC 3
49#define MIN_CONSTRAINED_BLEIC 4
50#define MIN_CONSTRAINED_BLIC 5
51#define MIN_CONSTRAINED_LEIC 6
52#define MIN_CONSTRAINED_LIC 7
53
54struct mexdata {
55 /* matlab names of the fitting function & its Jacobian */
56 char *fname, *jacname;
57
58 /* binary flags specifying if input p0 is a row or column vector */
59 int isrow_p0;
60
61 /* rhs args to be passed to matlab. rhs[0] is reserved for
62 * passing the parameter vector. If present, problem-specific
63 * data are passed in rhs[1], rhs[2], etc
64 */
65 mxArray **rhs;
66 int nrhs; /* >= 1 */
67};
68
69/* display printf-style error messages in matlab */
70static void matlabFmtdErrMsgTxt(char *fmt, ...)
71{
72char buf[256];
73va_list args;
74
75 va_start(args, fmt);
76 vsprintf(buf, fmt, args);
77 va_end(args);
78
79 mexErrMsgTxt(buf);
80}
81
82/* display printf-style warning messages in matlab */
83static void matlabFmtdWarnMsgTxt(char *fmt, ...)
84{
85char buf[256];
86va_list args;
87
88 va_start(args, fmt);
89 vsprintf(buf, fmt, args);
90 va_end(args);
91
92 mexWarnMsgTxt(buf);
93}
94
95static void func(double *p, double *hx, int m, int n, void *adata)
96{
97mxArray *lhs[1];
98double *mp, *mx;
99register int i;
100struct mexdata *dat=(struct mexdata *)adata;
101
102 /* prepare to call matlab */
103 mp=mxGetPr(dat->rhs[0]);
104 for(i=0; i<m; ++i)
105 mp[i]=p[i];
106
107 /* invoke matlab */
108 mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);
109
110 /* copy back results & cleanup */
111 mx=mxGetPr(lhs[0]);
112 for(i=0; i<n; ++i)
113 hx[i]=mx[i];
114
115 /* delete the matrix created by matlab */
116 mxDestroyArray(lhs[0]);
117}
118
119static void jacfunc(double *p, double *j, int m, int n, void *adata)
120{
121mxArray *lhs[1];
122double *mp;
123double *mj;
124register int i, k;
125struct mexdata *dat=(struct mexdata *)adata;
126
127 /* prepare to call matlab */
128 mp=mxGetPr(dat->rhs[0]);
129 for(i=0; i<m; ++i)
130 mp[i]=p[i];
131
132 /* invoke matlab */
133 mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
134
135 /* copy back results & cleanup. Note that the nxm Jacobian
136 * computed by matlab should be transposed so that
137 * levmar gets it in row major, as expected
138 */
139 mj=mxGetPr(lhs[0]);
140 for(i=0; i<n; ++i)
141 for(k=0; k<m; ++k)
142 j[i*m+k]=mj[i+k*n];
143
144 /* delete the matrix created by matlab */
145 mxDestroyArray(lhs[0]);
146}
147
148/* matlab matrices are in column-major, this routine converts them to row major for levmar */
149static double *getTranspose(mxArray *Am)
150{
151int m, n;
152double *At, *A;
153register int i, j;
154
155 m=mxGetM(Am);
156 n=mxGetN(Am);
157 A=mxGetPr(Am);
158 At=mxMalloc(m*n*sizeof(double));
159
160 for(i=0; i<m; i++)
161 for(j=0; j<n; j++)
162 At[i*n+j]=A[i+j*m];
163
164 return At;
165}
166
167/* check the supplied matlab function and its Jacobian. Returns 1 on error, 0 otherwise */
168static int checkFuncAndJacobian(double *p, int m, int n, int havejac, struct mexdata *dat)
169{
170mxArray *lhs[1];
171register int i;
172int ret=0;
173double *mp;
174
175 mexSetTrapFlag(1); /* handle errors in the MEX-file */
176
177 mp=mxGetPr(dat->rhs[0]);
178 for(i=0; i<m; ++i)
179 mp[i]=p[i];
180
181 /* attempt to call the supplied func */
182 i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);
183 if(i){
184 fprintf(stderr, "levmar: error calling '%s'.\n", dat->fname);
185 ret=1;
186 }
187 else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) ||
188 __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){
189 fprintf(stderr, "levmar: '%s' should produce a real vector with %d elements (got %d).\n",
190 dat->fname, n, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0])));
191 ret=1;
192 }
193 /* delete the matrix created by matlab */
194 mxDestroyArray(lhs[0]);
195
196 if(havejac){
197 /* attempt to call the supplied jac */
198 i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
199 if(i){
200 fprintf(stderr, "levmar: error calling '%s'.\n", dat->jacname);
201 ret=1;
202 }
203 else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || mxGetM(lhs[0])!=n || mxGetN(lhs[0])!=m){
204 fprintf(stderr, "levmar: '%s' should produce a real %dx%d matrix (got %dx%d).\n",
205 dat->jacname, n, m, mxGetM(lhs[0]), mxGetN(lhs[0]));
206 ret=1;
207 }
208 else if(mxIsSparse(lhs[0])){
209 fprintf(stderr, "levmar: '%s' should produce a real dense matrix (got a sparse one).\n", dat->jacname);
210 ret=1;
211 }
212 /* delete the matrix created by matlab */
213 mxDestroyArray(lhs[0]);
214 }
215
216 mexSetTrapFlag(0); /* on error terminate the MEX-file and return control to the MATLAB prompt */
217
218 return ret;
219}
220
221
222/*
223[ret, p, info, covar]=levmar_der (f, j, p0, x, itmax, opts, 'unc' ...)
224[ret, p, info, covar]=levmar_bc (f, j, p0, x, itmax, opts, 'bc', lb, ub, ...)
225[ret, p, info, covar]=levmar_bc (f, j, p0, x, itmax, opts, 'bc', lb, ub, dscl, ...)
226[ret, p, info, covar]=levmar_lec (f, j, p0, x, itmax, opts, 'lec', A, b, ...)
227[ret, p, info, covar]=levmar_blec(f, j, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...)
228
229[ret, p, info, covar]=levmar_bleic(f, j, p0, x, itmax, opts, 'bleic', lb, ub, A, b, C, d, ...)
230[ret, p, info, covar]=levmar_blic (f, j, p0, x, itmax, opts, 'blic', lb, ub, C, d, ...)
231[ret, p, info, covar]=levmar_leic (f, j, p0, x, itmax, opts, 'leic', A, b, C, d, ...)
232[ret, p, info, covar]=levmar_lic (f, j, p0, x, itmax, opts, 'lic', C, d, ...)
233
234*/
235
236void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *Prhs[])
237{
238register int i;
239register double *pdbl;
240mxArray **prhs=(mxArray **)&Prhs[0], *At, *Ct;
241struct mexdata mdata;
242int len, status;
243double *p, *p0, *ret, *x;
244int m, n, havejac, Arows, Crows, itmax, nopts, mintype, nextra;
245double opts[LM_OPTS_SZ]={LM_INIT_MU, LM_STOP_THRESH, LM_STOP_THRESH, LM_STOP_THRESH, LM_DIFF_DELTA};
246double info[LM_INFO_SZ];
247double *lb=NULL, *ub=NULL, *dscl=NULL, *A=NULL, *b=NULL, *wghts=NULL, *C=NULL, *d=NULL, *covar=NULL;
248
249 /* parse input args; start by checking their number */
250 if((nrhs<5))
251 matlabFmtdErrMsgTxt("levmar: at least 5 input arguments required (got %d).", nrhs);
252 if(nlhs>4)
253 matlabFmtdErrMsgTxt("levmar: too many output arguments (max. 4, got %d).", nlhs);
254 else if(nlhs<2)
255 matlabFmtdErrMsgTxt("levmar: too few output arguments (min. 2, got %d).", nlhs);
256
257 /* note that in order to accommodate optional args, prhs & nrhs are adjusted accordingly below */
258
259 /** func **/
260 /* first argument must be a string , i.e. a char row vector */
261 if(mxIsChar(prhs[0])!=1)
262 mexErrMsgTxt("levmar: first argument must be a string.");
263 if(mxGetM(prhs[0])!=1)
264 mexErrMsgTxt("levmar: first argument must be a string (i.e. char row vector).");
265 /* store supplied name */
266 len=mxGetN(prhs[0])+1;
267 mdata.fname=mxCalloc(len, sizeof(char));
268 status=mxGetString(prhs[0], mdata.fname, len);
269 if(status!=0)
270 mexErrMsgTxt("levmar: not enough space. String is truncated.");
271
272 /** jac (optional) **/
273 /* check whether second argument is a string */
274 if(mxIsChar(prhs[1])==1){
275 if(mxGetM(prhs[1])!=1)
276 mexErrMsgTxt("levmar: second argument must be a string (i.e. row vector).");
277 /* store supplied name */
278 len=mxGetN(prhs[1])+1;
279 mdata.jacname=mxCalloc(len, sizeof(char));
280 status=mxGetString(prhs[1], mdata.jacname, len);
281 if(status!=0)
282 mexErrMsgTxt("levmar: not enough space. String is truncated.");
283 havejac=1;
284
285 ++prhs;
286 --nrhs;
287 }
288 else{
289 mdata.jacname=NULL;
290 havejac=0;
291 }
292
293#ifdef DEBUG
294 fflush(stderr);
295 fprintf(stderr, "LEVMAR: %s analytic Jacobian\n", havejac? "with" : "no");
296#endif /* DEBUG */
297
298/* CHECK
299if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGetN(prhs[1])==1))
300*/
301
302 /** p0 **/
303 /* the second required argument must be a real row or column vector */
304 if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 || mxGetN(prhs[1])==1))
305 mexErrMsgTxt("levmar: p0 must be a real vector.");
306 p0=mxGetPr(prhs[1]);
307 /* determine if we have a row or column vector and retrieve its
308 * size, i.e. the number of parameters
309 */
310 if(mxGetM(prhs[1])==1){
311 m=mxGetN(prhs[1]);
312 mdata.isrow_p0=1;
313 }
314 else{
315 m=mxGetM(prhs[1]);
316 mdata.isrow_p0=0;
317 }
318 /* copy input parameter vector to avoid destroying it */
319 p=mxMalloc(m*sizeof(double));
320 for(i=0; i<m; ++i)
321 p[i]=p0[i];
322
323 /** x **/
324 /* the third required argument must be a real row or column vector */
325 if(!mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]) || !(mxGetM(prhs[2])==1 || mxGetN(prhs[2])==1))
326 mexErrMsgTxt("levmar: x must be a real vector.");
327 x=mxGetPr(prhs[2]);
328 n=__MAX__(mxGetM(prhs[2]), mxGetN(prhs[2]));
329
330 /** itmax **/
331 /* the fourth required argument must be a scalar */
332 if(!mxIsDouble(prhs[3]) || mxIsComplex(prhs[3]) || mxGetM(prhs[3])!=1 || mxGetN(prhs[3])!=1)
333 mexErrMsgTxt("levmar: itmax must be a scalar.");
334 itmax=(int)mxGetScalar(prhs[3]);
335
336 /** opts **/
337 /* the fifth required argument must be a real row or column vector */
338 if(!mxIsDouble(prhs[4]) || mxIsComplex(prhs[4]) || (!(mxGetM(prhs[4])==1 || mxGetN(prhs[4])==1) &&
339 !(mxGetM(prhs[4])==0 && mxGetN(prhs[4])==0)))
340 mexErrMsgTxt("levmar: opts must be a real vector.");
341 pdbl=mxGetPr(prhs[4]);
342 nopts=__MAX__(mxGetM(prhs[4]), mxGetN(prhs[4]));
343 if(nopts!=0){ /* if opts==[], nothing needs to be done and the defaults are used */
344 if(nopts>LM_OPTS_SZ)
345 matlabFmtdErrMsgTxt("levmar: opts must have at most %d elements, got %d.", LM_OPTS_SZ, nopts);
346 else if(nopts<((havejac)? LM_OPTS_SZ-1 : LM_OPTS_SZ))
347 matlabFmtdWarnMsgTxt("levmar: only the %d first elements of opts specified, remaining set to defaults.", nopts);
348 for(i=0; i<nopts; ++i)
349 opts[i]=pdbl[i];
350 }
351#ifdef DEBUG
352 else{
353 fflush(stderr);
354 fprintf(stderr, "LEVMAR: empty options vector, using defaults\n");
355 }
356#endif /* DEBUG */
357
358 /** mintype (optional) **/
359 /* check whether sixth argument is a string */
360 if(nrhs>=6 && mxIsChar(prhs[5])==1 && mxGetM(prhs[5])==1){
361 char *minhowto;
362
363 /* examine supplied name */
364 len=mxGetN(prhs[5])+1;
365 minhowto=mxCalloc(len, sizeof(char));
366 status=mxGetString(prhs[5], minhowto, len);
367 if(status!=0)
368 mexErrMsgTxt("levmar: not enough space. String is truncated.");
369
370 for(i=0; minhowto[i]; ++i)
371 minhowto[i]=tolower(minhowto[i]);
372 if(!strncmp(minhowto, "unc", 3)) mintype=MIN_UNCONSTRAINED;
373 else if(!strncmp(minhowto, "bc", 2)) mintype=MIN_CONSTRAINED_BC;
374 else if(!strncmp(minhowto, "lec", 3)) mintype=MIN_CONSTRAINED_LEC;
375 else if(!strncmp(minhowto, "blec", 4)) mintype=MIN_CONSTRAINED_BLEC;
376 else if(!strncmp(minhowto, "bleic", 5)) mintype=MIN_CONSTRAINED_BLEIC;
377 else if(!strncmp(minhowto, "blic", 4)) mintype=MIN_CONSTRAINED_BLIC;
378 else if(!strncmp(minhowto, "leic", 4)) mintype=MIN_CONSTRAINED_LEIC;
379 else if(!strncmp(minhowto, "lic", 3)) mintype=MIN_CONSTRAINED_BLIC;
380 else matlabFmtdErrMsgTxt("levmar: unknown minimization type '%s'.", minhowto);
381
382 mxFree(minhowto);
383
384 ++prhs;
385 --nrhs;
386 }
387 else
388 mintype=MIN_UNCONSTRAINED;
389
390 if(mintype==MIN_UNCONSTRAINED) goto extraargs;
391
392 /* arguments below this point are optional and their presence depends
393 * upon the minimization type determined above
394 */
395 /** lb, ub **/
396 if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_BLEIC)){
397 /* check if the next two arguments are real row or column vectors */
398 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
399 if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
400 if((i=__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5])))!=m)
401 matlabFmtdErrMsgTxt("levmar: lb must have %d elements, got %d.", m, i);
402 if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=m)
403 matlabFmtdErrMsgTxt("levmar: ub must have %d elements, got %d.", m, i);
404
405 lb=mxGetPr(prhs[5]);
406 ub=mxGetPr(prhs[6]);
407
408 prhs+=2;
409 nrhs-=2;
410 }
411 }
412 }
413
414 /** dscl **/
415 if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC)){
416 /* check if the next argument is a real row or column vector */
417 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
418 if((i=__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5])))!=m)
419 matlabFmtdErrMsgTxt("levmar: dscl must have %d elements, got %d.", m, i);
420
421 dscl=mxGetPr(prhs[5]);
422
423 ++prhs;
424 --nrhs;
425 }
426 }
427
428 /** A, b **/
429 if(nrhs>=7 && (mintype==MIN_CONSTRAINED_LEC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_BLEIC)){
430 /* check if the next two arguments are a real matrix and a real row or column vector */
431 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){
432 if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
433 if((i=mxGetN(prhs[5]))!=m)
434 matlabFmtdErrMsgTxt("levmar: A must have %d columns, got %d.", m, i);
435 if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=(Arows=mxGetM(prhs[5])))
436 matlabFmtdErrMsgTxt("levmar: b must have %d elements, got %d.", Arows, i);
437
438 At=prhs[5];
439 b=mxGetPr(prhs[6]);
440 A=getTranspose(At);
441
442 prhs+=2;
443 nrhs-=2;
444 }
445 }
446 }
447
448 /* wghts */
449 /* check if we have a weights vector */
450 if(nrhs>=6 && mintype==MIN_CONSTRAINED_BLEC){ /* only check if we have seen both box & linear constraints */
451 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
452 if(__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5]))==m){
453 wghts=mxGetPr(prhs[5]);
454
455 ++prhs;
456 --nrhs;
457 }
458 }
459 }
460
461 /** C, d **/
462 if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BLEIC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_LIC)){
463 /* check if the next two arguments are a real matrix and a real row or column vector */
464 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){
465 if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
466 if((i=mxGetN(prhs[5]))!=m)
467 matlabFmtdErrMsgTxt("levmar: C must have %d columns, got %d.", m, i);
468 if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=(Crows=mxGetM(prhs[5])))
469 matlabFmtdErrMsgTxt("levmar: d must have %d elements, got %d.", Crows, i);
470
471 Ct=prhs[5];
472 d=mxGetPr(prhs[6]);
473 C=getTranspose(Ct);
474
475 prhs+=2;
476 nrhs-=2;
477 }
478 }
479 }
480
481 /* arguments below this point are assumed to be extra arguments passed
482 * to every invocation of the fitting function and its Jacobian
483 */
484
485extraargs:
486 /* handle any extra args and allocate memory for
487 * passing the current parameter estimate to matlab
488 */
489 nextra=nrhs-5;
490 mdata.nrhs=nextra+1;
491 mdata.rhs=(mxArray **)mxMalloc(mdata.nrhs*sizeof(mxArray *));
492 for(i=0; i<nextra; ++i)
493 mdata.rhs[i+1]=(mxArray *)prhs[nrhs-nextra+i]; /* discard 'const' modifier */
494#ifdef DEBUG
495 fflush(stderr);
496 fprintf(stderr, "LEVMAR: %d extra args\n", nextra);
497#endif /* DEBUG */
498
499 if(mdata.isrow_p0){ /* row vector */
500 mdata.rhs[0]=mxCreateDoubleMatrix(1, m, mxREAL);
501 /*
502 mxSetM(mdata.rhs[0], 1);
503 mxSetN(mdata.rhs[0], m);
504 */
505 }
506 else{ /* column vector */
507 mdata.rhs[0]=mxCreateDoubleMatrix(m, 1, mxREAL);
508 /*
509 mxSetM(mdata.rhs[0], m);
510 mxSetN(mdata.rhs[0], 1);
511 */
512 }
513
514 /* ensure that the supplied function & Jacobian are as expected */
515 if(checkFuncAndJacobian(p, m, n, havejac, &mdata)){
516 status=LM_ERROR;
517 goto cleanup;
518 }
519
520 if(nlhs>3) /* covariance output required */
521 covar=mxMalloc(m*m*sizeof(double));
522
523 /* invoke levmar */
524 switch(mintype){
525 case MIN_UNCONSTRAINED: /* no constraints */
526 if(havejac)
527 status=dlevmar_der(func, jacfunc, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
528 else
529 status=dlevmar_dif(func, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
530#ifdef DEBUG
531 fflush(stderr);
532 fprintf(stderr, "LEVMAR: calling dlevmar_der()/dlevmar_dif()\n");
533#endif /* DEBUG */
534 break;
535 case MIN_CONSTRAINED_BC: /* box constraints */
536 if(havejac)
537 status=dlevmar_bc_der(func, jacfunc, p, x, m, n, lb, ub, dscl, itmax, opts, info, NULL, covar, (void *)&mdata);
538 else
539 status=dlevmar_bc_dif(func, p, x, m, n, lb, ub, dscl, itmax, opts, info, NULL, covar, (void *)&mdata);
540#ifdef DEBUG
541 fflush(stderr);
542 fprintf(stderr, "LEVMAR: calling dlevmar_bc_der()/dlevmar_bc_dif()\n");
543#endif /* DEBUG */
544 break;
545 case MIN_CONSTRAINED_LEC: /* linear equation constraints */
546#ifdef HAVE_LAPACK
547 if(havejac)
548 status=dlevmar_lec_der(func, jacfunc, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata);
549 else
550 status=dlevmar_lec_dif(func, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata);
551#else
552 mexErrMsgTxt("levmar: no linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
553#endif /* HAVE_LAPACK */
554
555#ifdef DEBUG
556 fflush(stderr);
557 fprintf(stderr, "LEVMAR: calling dlevmar_lec_der()/dlevmar_lec_dif()\n");
558#endif /* DEBUG */
559 break;
560 case MIN_CONSTRAINED_BLEC: /* box & linear equation constraints */
561#ifdef HAVE_LAPACK
562 if(havejac)
563 status=dlevmar_blec_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
564 else
565 status=dlevmar_blec_dif(func, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
566#else
567 mexErrMsgTxt("levmar: no box & linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
568#endif /* HAVE_LAPACK */
569
570#ifdef DEBUG
571 fflush(stderr);
572 fprintf(stderr, "LEVMAR: calling dlevmar_blec_der()/dlevmar_blec_dif()\n");
573#endif /* DEBUG */
574 break;
575 case MIN_CONSTRAINED_BLEIC: /* box, linear equation & inequalities constraints */
576#ifdef HAVE_LAPACK
577 if(havejac)
578 status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
579 else
580 status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
581#else
582 mexErrMsgTxt("levmar: no box, linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
583#endif /* HAVE_LAPACK */
584
585#ifdef DEBUG
586 fflush(stderr);
587 fprintf(stderr, "LEVMAR: calling dlevmar_bleic_der()/dlevmar_bleic_dif()\n");
588#endif /* DEBUG */
589 break;
590 case MIN_CONSTRAINED_BLIC: /* box, linear inequalities constraints */
591#ifdef HAVE_LAPACK
592 if(havejac)
593 status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
594 else
595 status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
596#else
597 mexErrMsgTxt("levmar: no box & linear inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
598#endif /* HAVE_LAPACK */
599
600#ifdef DEBUG
601 fflush(stderr);
602 fprintf(stderr, "LEVMAR: calling dlevmar_blic_der()/dlevmar_blic_dif()\n");
603#endif /* DEBUG */
604 break;
605 case MIN_CONSTRAINED_LEIC: /* linear equation & inequalities constraints */
606#ifdef HAVE_LAPACK
607 if(havejac)
608 status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
609 else
610 status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
611#else
612 mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
613#endif /* HAVE_LAPACK */
614
615#ifdef DEBUG
616 fflush(stderr);
617 fprintf(stderr, "LEVMAR: calling dlevmar_leic_der()/dlevmar_leic_dif()\n");
618#endif /* DEBUG */
619 break;
620 case MIN_CONSTRAINED_LIC: /* linear inequalities constraints */
621#ifdef HAVE_LAPACK
622 if(havejac)
623 status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
624 else
625 status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata);
626#else
627 mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
628#endif /* HAVE_LAPACK */
629
630#ifdef DEBUG
631 fflush(stderr);
632 fprintf(stderr, "LEVMAR: calling dlevmar_lic_der()/dlevmar_lic_dif()\n");
633#endif /* DEBUG */
634 break;
635 default:
636 mexErrMsgTxt("levmar: unexpected internal error.");
637 }
638
639#ifdef DEBUG
640 fflush(stderr);
641 printf("LEVMAR: minimization returned %d in %g iter, reason %g\n\tSolution: ", status, info[5], info[6]);
642 for(i=0; i<m; ++i)
643 printf("%.7g ", p[i]);
644 printf("\n\n\tMinimization info:\n\t");
645 for(i=0; i<LM_INFO_SZ; ++i)
646 printf("%g ", info[i]);
647 printf("\n");
648#endif /* DEBUG */
649
650 /* copy back return results */
651 /** ret **/
652 plhs[0]=mxCreateDoubleMatrix(1, 1, mxREAL);
653 ret=mxGetPr(plhs[0]);
654 ret[0]=(double)status;
655
656 /** popt **/
657 plhs[1]=(mdata.isrow_p0==1)? mxCreateDoubleMatrix(1, m, mxREAL) : mxCreateDoubleMatrix(m, 1, mxREAL);
658 pdbl=mxGetPr(plhs[1]);
659 for(i=0; i<m; ++i)
660 pdbl[i]=p[i];
661
662 /** info **/
663 if(nlhs>2){
664 plhs[2]=mxCreateDoubleMatrix(1, LM_INFO_SZ, mxREAL);
665 pdbl=mxGetPr(plhs[2]);
666 for(i=0; i<LM_INFO_SZ; ++i)
667 pdbl[i]=info[i];
668 }
669
670 /** covar **/
671 if(nlhs>3){
672 plhs[3]=mxCreateDoubleMatrix(m, m, mxREAL);
673 pdbl=mxGetPr(plhs[3]);
674 for(i=0; i<m*m; ++i) /* covariance matrices are symmetric, thus no need to transpose! */
675 pdbl[i]=covar[i];
676 }
677
678cleanup:
679 /* cleanup */
680 mxDestroyArray(mdata.rhs[0]);
681 if(A) mxFree(A);
682 if(C) mxFree(C);
683
684 mxFree(mdata.fname);
685 if(havejac) mxFree(mdata.jacname);
686 mxFree(p);
687 mxFree(mdata.rhs);
688 if(covar) mxFree(covar);
689
690 if(status==LM_ERROR)
691 mexWarnMsgTxt("levmar: optimization returned with an error!");
692}
Note: See TracBrowser for help on using the repository browser.