Line data Source code
1 : /** \file levmarInterface.hpp
2 : * \author Jared R. Males
3 : * \brief A c++ interface to the templatized levmar minimization routines..
4 : * \ingroup fitting_files
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2015, 2016, 2017 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef levmarInterface_hpp
28 : #define levmarInterface_hpp
29 :
30 : #include <iostream>
31 :
32 : #include "../../mxlib.hpp"
33 :
34 : #include "templateLevmar.hpp"
35 : #include "../../sys/timeUtils.hpp"
36 :
37 : namespace mx
38 : {
39 : namespace math
40 : {
41 : namespace fit
42 : {
43 :
44 : // Forwards
45 : template <typename T>
46 : struct hasJacobian;
47 :
48 : /// A templatized interface to the levmar package
49 : /** Requires a fitter class, which conforms to one of the following minimum specifications.
50 : * To use the finite difference jacobian calculation:
51 : * \code
52 : * //This will cause the levmar_dif routine to be used
53 : * template<typename _realT>
54 : * struct dif_fitter
55 : * {
56 : * typedef _realT realT; //required
57 : *
58 : * static void func( realT *p, // [in] the parameter values, size m
59 : * realT *hx, // [out] the error for this p, size n
60 : * int m, // [in] the number of parameters
61 : * int n, // [in] the number of data points
62 : * void *adata // [in] auxiliary data
63 : * )
64 : * {
65 : * //do stuff here . . .
66 : * }
67 : * };
68 : * \endcode
69 : * If you wish to provide your own jacobian:
70 : * \code
71 : * //This will cause the levmar_der routine to be used
72 : * template<typename _realT>
73 : * struct der_fitter
74 : * {
75 : * typedef _realT realT; //required
76 : *
77 : * typdef bool hasJacobian; //this signals that jacf exists and should be used.
78 : *
79 : * static void func( realT *p, // [in] the parameter values, size m
80 : * realT *hx, // [out] the error for this p, size n
81 : * int m, // [in] the number of parameters
82 : * int n, // [in] the number of data points
83 : * void *adata // [in] auxiliary data
84 : * )
85 : * {
86 : * //do stuff here . . .
87 : * }
88 : *
89 : * // This is the jacobian.
90 : * static void jacf( realT *p, // [in] the parameter values, size m
91 : * realT *j, // [out] the jacobian for this p, size n
92 : * int m, // [in] the number of parameters
93 : * int n, // [in] the number of data points
94 : * void *adata // [in] auxiliary data
95 : * )
96 : * {
97 : * //do stuff here . . .
98 : * }
99 : *
100 : * };
101 : *
102 : * \endcode
103 : *
104 : * Note that if your fitter has a jacobian function (jacf), you must define
105 : * \code
106 : * typedef bool hasJacobian;
107 : * \endcode
108 : * for it to be used.
109 : *
110 : * \tparam fitterT a class with at least the minimum interface described above.
111 : *
112 : * \ingroup fitting
113 : */
114 : template <class fitterT>
115 : class levmarInterface
116 : {
117 :
118 : public:
119 : typedef typename fitterT::realT realT;
120 :
121 : protected:
122 : realT *p; ///< Parameter array. On input is the initial estimates. On output has the estimated solution.
123 : realT *init_p; ///< Parameter array on input, saved for comparison.
124 :
125 : int m; ///< Parameter vector dimension (i.e. number of unknowns)
126 : bool own_p; ///< Flag indicating whether the p array is owned by this object (for de-allocation).
127 :
128 : int itmax; ///< Maximum number of iterations, default is 100
129 :
130 : public:
131 : realT *x; ///< I: measurement vector. NULL implies a zero vector
132 :
133 : int n; ///< I: measurement vector dimension
134 :
135 : /// Options passed to the minimization routines. See \ref set_opts for details.
136 : realT opts[LM_OPTS_SZ];
137 :
138 : /// Information regarding the minimization.
139 : /** See the levmar source code for documentation. These fields are accessed by
140 : * \ref get_initial_norm, \ref get_final_norm, \ref get_iterations, \ref get_reason_code,
141 : * \ref get_reason_string, \ref get_fevals, \ref get_jevals, \ref get_nlinsys
142 : */
143 : realT info[LM_INFO_SZ];
144 :
145 : /// Elapsed time of the fitting procedure
146 : double deltaT;
147 :
148 : /// Working memory passed to the levmar routines.
149 : /** From the levmar documentation: at least LM_DER/DIF_WORKSZ() reals large, allocated if NULL
150 : * Here this is always allocated by a call to \ref allocate_work.
151 : */
152 : realT *work;
153 :
154 : /// The current size of the work array
155 : int work_sz;
156 :
157 : /// Covariance matrix corresponding to LS solution; mxm.
158 : /** Here this is allocated to size m-x-m by allocate_work, but if you don't want it
159 : * allocated and calculated
160 : * set \ref getCovar to false and NULL will be passed.
161 : */
162 : realT *covar;
163 :
164 : /// The current size of the covar array
165 : int covar_sz;
166 :
167 : /// Controls whether the covar array is allocated.
168 : bool getCovar;
169 :
170 : /// Pointer to possibly additional data, passed uninterpreted to func & jacf.
171 : /** Set to NULL if not needed.
172 : */
173 : void *adata;
174 :
175 : private:
176 : /// Initialization common to all constructors
177 : void initialize();
178 :
179 : public:
180 : /// Default constructor.
181 : /** With this constructor, you must set the parameter, data, and adata before calling fit().
182 : */
183 : levmarInterface();
184 :
185 : /// Setup constructor
186 : /** with this constructor fit() will work immediately.
187 : *
188 : */
189 : levmarInterface( realT *i_p, ///< [in] pointer to the initial parameter guess
190 : realT *i_x, ///< [in] pointer to the data (can be NULL)
191 : int i_m, ///< [in] the size of i_p
192 : int i_n, ///< [in] the size of i_x
193 : void *i_adata ///< [in] pointer to auxiliary data (can be NULL)
194 : );
195 :
196 : /// Destructor
197 : /** Frees the work and covar matrices.
198 : */
199 : ~levmarInterface();
200 :
201 : /// Set number of parameters, but don't allocate
202 : void nParams( int i_m /**< [in] the number of parameters */ );
203 :
204 : /// Get the current number of parameters
205 : /** \returns the current number of parameters (m)
206 : */
207 : int nParams();
208 :
209 : /// Allocate parameters array based on previous call to \ref nParams
210 : void allocate_params();
211 :
212 : /// Set number of parameters and allocate
213 : void allocate_params( int i_m );
214 :
215 : /// Point the parameter pointer at an externally allocated array
216 : void point_params( realT *i_p );
217 :
218 : /// Point the parameter pointer at an externally allocated array
219 : void point_params( realT *i_p, int i_m );
220 :
221 : /// Copy parameters to the parameter array
222 : /** This assumes that either the array was allocated (i.e. with \ref allocate_params) or
223 : * that \ref point_params has been called
224 : */
225 : void set_params( realT *i_p );
226 :
227 : /// Get current pointer array address
228 : realT *get_params();
229 :
230 : /// Set the maximum number of iterations
231 : /** Sets itmax. Initialization default is itmax = 100
232 : *
233 : * \param i_itmax the new value of itmax to set
234 : */
235 : void set_itmax( int i_itmax );
236 :
237 : /// Get the maximum number of iterations
238 : int get_itmax();
239 :
240 : /// Allocate the work and covar matrices
241 : /** Uses a function object specialized for whether or not there is a jacobian
242 : * to determine the size of work.
243 : */
244 : void allocate_work();
245 :
246 : /// Set one of the minimization options to val
247 : /** The options correspond to:
248 : *
249 : * 0: the scale factor of the initial \f$ \mu \f$, default is 1e-3.
250 : *
251 : * 1: \f$ \epsilon_1 \f$ Stopping threshold for ||J^T e||_inf, default is 1e-17.
252 : *
253 : * 2: \f$ \epsilon_2 \f$ Stopping threshold for ||Dp||_2, default is 1e-17.
254 : *
255 : * 3: \f$ \epsilon_3 \f$ Stopping threshold for ||e||_2, default is 1e-17.
256 : *
257 : * 4: \f$ \Delta_{diff}\f$ Stepsize for finite differences, default is 1e-6.
258 : *
259 : * \param n is the option number
260 : * \param val is the value to set
261 : */
262 : void set_opts( int n, realT val );
263 :
264 : /// Set one or all of the minimization options to the default.
265 : /** See \ref set_opts for discription of the options.
266 : *
267 : * \param n the option number. Pass -1 to set all options to defaults.
268 : */
269 : void set_opts_default( int n = -1 );
270 :
271 : /// Get the current value of an option.
272 : /** See \ref set_opts for a description of the options
273 : *
274 : * \param n the option number
275 : */
276 : realT get_opts( int n );
277 :
278 : /// Perform the fit
279 : /** This calls \ref allocate_work, and then dispatches the levmar routine appropriate for fitterT.
280 : */
281 : int fit();
282 :
283 : /// Returns the L2-norm before minimization occurs
284 : realT get_initial_norm();
285 :
286 : /// Returns the L2-norm at the end of the minimization
287 : realT get_final_norm();
288 :
289 : /// Get the number of iterations taken during the minimization
290 : int get_iterations();
291 :
292 : /// Get a code specifying the reason minimization terminated.
293 : int get_reason_code();
294 :
295 : /// Get the descriptive string describing the reason minimization terminated.
296 : std::string get_reason_string();
297 :
298 : /// Get the number of function evaluations during the minimization
299 : int get_fevals();
300 :
301 : /// Get the number of jacobian evaluations during the minimization
302 : int get_jevals();
303 :
304 : /// Get the number of linear system solutions during the minimization
305 : int get_nlinsys();
306 :
307 : /// Get the elapsed time of the fit
308 : double get_deltaT()
309 : {
310 : return deltaT;
311 : }
312 :
313 : // Status reports
314 :
315 : /// Output current parameters to a stream
316 : /** Prints a formatted list of all current fit parameters.
317 : *
318 : * \tparam iosT is a std::ostream-like type.
319 : * \tparam comment is a comment character to start each line. Can be '\0'.
320 : */
321 : template <typename iosT, char comment = '#'>
322 : iosT &dumpParameters( iosT &ios /**< [in] a std::ostream-like stream. */ );
323 :
324 : /// Dump the parameter vector to stdout.
325 : /**
326 : * \tparam comment is a comment character to start each line. Can be '\0'.
327 : */
328 : template <char comment = '#'>
329 : std::ostream &dumpParameters();
330 :
331 : /// Output current parameters to a stream
332 : /** Prints a formatted list of all current fit parameters.
333 : *
334 : * \tparam iosT is a std::ostream-like type.
335 : * \tparam comment is a comment character to start each line. Can be '\0'.
336 : */
337 : template <typename iosT, char comment = '#'>
338 : iosT &dumpReport( iosT &ios, ///< [in] a std::ostream-like stream.
339 : bool dumpParams = true ///< [in] [optional] whether or not to dump the parameters.
340 : );
341 :
342 : /// Dump a status report to stdout
343 : /**
344 : * \tparam comment is a comment character to start each line. Can be '\0'.
345 : */
346 : template <char comment = '#'>
347 : std::ostream &dumpReport( bool dumpParams = true /**< [in] [optional] whether or not to dump the parameters.*/ );
348 : };
349 :
350 : template <class fitterT>
351 0 : void levmarInterface<fitterT>::initialize()
352 : {
353 0 : p = 0;
354 0 : own_p = false;
355 :
356 0 : init_p = 0;
357 :
358 0 : m = 0;
359 :
360 0 : x = 0;
361 0 : n = 0;
362 0 : set_opts_default( -1 ); // set all opts to defaults.
363 0 : work = 0;
364 0 : work_sz = 0;
365 0 : covar = 0;
366 0 : covar_sz = 0;
367 0 : getCovar = true;
368 0 : adata = 0;
369 :
370 0 : for( int i = 0; i < LM_INFO_SZ; ++i )
371 0 : info[i] = 0;
372 0 : deltaT = 0;
373 :
374 0 : itmax = 100;
375 0 : }
376 :
377 : template <class fitterT>
378 0 : levmarInterface<fitterT>::levmarInterface()
379 : {
380 0 : initialize();
381 0 : }
382 :
383 : template <class fitterT>
384 : levmarInterface<fitterT>::levmarInterface(
385 : typename fitterT::realT *i_p, typename fitterT::realT *i_x, int i_m, int i_n, void *i_adata )
386 : {
387 : initialize();
388 :
389 : p = i_p;
390 : x = i_x;
391 : m = i_m;
392 : n = i_n;
393 : adata = i_adata;
394 : }
395 :
396 : template <class fitterT>
397 0 : levmarInterface<fitterT>::~levmarInterface()
398 : {
399 0 : if( p && own_p )
400 0 : free( p );
401 :
402 0 : if( init_p )
403 0 : free( init_p );
404 :
405 0 : if( work )
406 0 : free( work );
407 :
408 0 : if( covar )
409 0 : free( covar );
410 0 : }
411 :
412 : template <class fitterT>
413 : void levmarInterface<fitterT>::nParams( int i_m )
414 : {
415 : // If we own and have allocated p, then de-alloc
416 : if( p && own_p )
417 : {
418 : free( p );
419 : p = 0;
420 : }
421 :
422 : m = i_m;
423 :
424 : // Also allocate the init_p storage for initial guess
425 : if( init_p )
426 : {
427 : free( init_p );
428 : }
429 :
430 : init_p = (typename fitterT::realT *)malloc( sizeof( typename fitterT::realT ) * m );
431 : }
432 :
433 : template <class fitterT>
434 : int levmarInterface<fitterT>::nParams()
435 : {
436 : return m;
437 : }
438 :
439 : template <class fitterT>
440 : void levmarInterface<fitterT>::allocate_params()
441 : {
442 : if( p && own_p )
443 : {
444 : free( p );
445 : }
446 :
447 : p = (typename fitterT::realT *)malloc( sizeof( typename fitterT::realT ) * m );
448 :
449 : own_p = true;
450 : }
451 :
452 : template <class fitterT>
453 : void levmarInterface<fitterT>::allocate_params( int i_m )
454 : {
455 : nParams( i_m );
456 :
457 : allocate_params();
458 : }
459 :
460 : template <class fitterT>
461 : void levmarInterface<fitterT>::point_params( realT *i_p )
462 : {
463 : if( p && own_p )
464 : {
465 : free( p );
466 : }
467 :
468 : p = i_p;
469 : }
470 :
471 : template <class fitterT>
472 : void levmarInterface<fitterT>::point_params( realT *i_p, int i_m )
473 : {
474 : nParams( i_m );
475 : point_params( i_p );
476 : }
477 :
478 : template <class fitterT>
479 : void levmarInterface<fitterT>::set_itmax( int i_itmax )
480 : {
481 : itmax = i_itmax;
482 : }
483 :
484 : template <class fitterT>
485 : int levmarInterface<fitterT>::get_itmax()
486 : {
487 : return itmax;
488 : }
489 :
490 : template <class fitterT>
491 : typename fitterT::realT *levmarInterface<fitterT>::get_params()
492 : {
493 : return p;
494 : }
495 :
496 : template <class fitterT>
497 : void levmarInterface<fitterT>::set_params( realT *i_p )
498 : {
499 : for( int i = 0; i < m; i++ )
500 : p[i] = i_p[i];
501 : }
502 :
503 : // Functor which is used by allocate() if hasJacobian is false
504 : template <class fitterT, bool jacf = hasJacobian<fitterT>::value>
505 : struct levmar_allocate_size
506 : {
507 : typedef typename fitterT::realT realT;
508 :
509 : int operator()( int m, int n )
510 : {
511 : return LM_DIF_WORKSZ( m, n );
512 : }
513 : };
514 :
515 : // Functor which is used by allocate() if hasJacobian is true
516 : template <class fitterT>
517 : struct levmar_allocate_size<fitterT, true>
518 : {
519 : typedef typename fitterT::realT realT;
520 :
521 : int operator()( int m, int n )
522 : {
523 : return LM_DER_WORKSZ( m, n );
524 : }
525 : };
526 :
527 : template <class fitterT>
528 : void levmarInterface<fitterT>::allocate_work()
529 : {
530 : // Create function object to get allocation size, which depends on whether there is Jacobian.
531 : levmar_allocate_size<fitterT> alloc_sz;
532 :
533 : if( work_sz < alloc_sz( m, n ) || !work )
534 : {
535 : if( work )
536 : free( work );
537 :
538 : work_sz = alloc_sz( m, n );
539 : work = (realT *)malloc( work_sz * sizeof( realT ) );
540 : }
541 :
542 : // Allocate if covar is desired and unallocated.
543 : if( getCovar )
544 : {
545 : if( covar_sz < m * m || !covar )
546 : {
547 : if( covar )
548 : free( covar );
549 :
550 : covar_sz = m * m;
551 : covar = (realT *)malloc( covar_sz * sizeof( realT ) );
552 : }
553 : }
554 : else
555 : {
556 : // If covar is not desired, de-allocate if allocated before.
557 : if( covar )
558 : free( covar );
559 : covar = 0;
560 : }
561 : }
562 :
563 : template <class fitterT>
564 : void levmarInterface<fitterT>::set_opts( int n, realT val )
565 : {
566 : opts[n] = val;
567 : }
568 :
569 : template <class fitterT>
570 0 : void levmarInterface<fitterT>::set_opts_default( int n )
571 : {
572 : // See the source in lm_core.c
573 :
574 0 : if( n == 0 || n == -1 )
575 : {
576 0 : opts[0] = LM_INIT_MU;
577 : }
578 :
579 0 : if( n > 0 && n < 4 )
580 : {
581 0 : opts[n] = LM_STOP_THRESH;
582 : }
583 0 : else if( n == -1 )
584 : {
585 0 : opts[1] = LM_STOP_THRESH;
586 0 : opts[2] = LM_STOP_THRESH;
587 0 : opts[3] = LM_STOP_THRESH;
588 : }
589 :
590 0 : if( n == 4 || n == -1 )
591 : {
592 0 : opts[4] = LM_DIFF_DELTA;
593 : }
594 0 : }
595 :
596 : template <class fitterT>
597 : typename fitterT::realT levmarInterface<fitterT>::get_opts( int n )
598 : {
599 : return opts[n];
600 : }
601 :
602 : // Functor which is used by fit() if hasJacobian is false
603 : template <class fitterT, bool jacf = hasJacobian<fitterT>::value>
604 : struct do_levmar
605 : {
606 : typedef typename fitterT::realT realT;
607 :
608 : int operator()(
609 : realT *p, realT *x, int m, int n, int itmax, realT *opts, realT *info, realT *work, realT *covar, void *adata )
610 : {
611 : return levmar_dif<realT>( &fitterT::func, p, x, m, n, itmax, opts, info, work, covar, adata );
612 : }
613 : };
614 :
615 : // Functor which is used by fit() if hasJacobian is true
616 : template <class fitterT>
617 : struct do_levmar<fitterT, true>
618 : {
619 : typedef typename fitterT::realT realT;
620 :
621 : int operator()(
622 : realT *p, realT *x, int m, int n, int itmax, realT *opts, realT *info, realT *work, realT *covar, void *adata )
623 : {
624 : return levmar_der<realT>( &fitterT::func, &fitterT::jacf, p, x, m, n, itmax, opts, info, work, covar, adata );
625 : }
626 : };
627 :
628 : template <class fitterT>
629 : int levmarInterface<fitterT>::fit()
630 : {
631 : realT *_opts;
632 :
633 : allocate_work();
634 :
635 : if( opts[0] == 0 )
636 : _opts = 0;
637 : else
638 : _opts = opts;
639 :
640 : // These may not be updated by the levmar library
641 : info[8] = 0;
642 : info[9] = 0;
643 :
644 : if( !init_p )
645 : {
646 : init_p = (typename fitterT::realT *)malloc( sizeof( typename fitterT::realT ) * m );
647 : }
648 :
649 : for( int i = 0; i < m; ++i )
650 : {
651 : init_p[i] = p[i];
652 : }
653 :
654 : double t0 = sys::get_curr_time();
655 :
656 : // Create one of the above functors, which depends on whether fitterT has a Jacobian.
657 : do_levmar<fitterT> fitter;
658 :
659 : fitter( p, x, m, n, itmax, _opts, info, work, covar, adata );
660 :
661 : deltaT = sys::get_curr_time() - t0;
662 :
663 : return 0;
664 : }
665 :
666 : template <class fitterT>
667 : typename fitterT::realT levmarInterface<fitterT>::get_initial_norm()
668 : {
669 : return info[0];
670 : }
671 :
672 : template <class fitterT>
673 : typename fitterT::realT levmarInterface<fitterT>::get_final_norm()
674 : {
675 : return info[1];
676 : }
677 :
678 : template <class fitterT>
679 : int levmarInterface<fitterT>::get_iterations()
680 : {
681 : return (int)info[5];
682 : }
683 :
684 : template <class fitterT>
685 : int levmarInterface<fitterT>::get_reason_code()
686 : {
687 : return (int)info[6];
688 : }
689 :
690 : template <class fitterT>
691 : int levmarInterface<fitterT>::get_fevals()
692 : {
693 : return (int)info[7];
694 : }
695 :
696 : template <class fitterT>
697 : int levmarInterface<fitterT>::get_jevals()
698 : {
699 : return (int)info[8];
700 : }
701 :
702 : template <class fitterT>
703 : int levmarInterface<fitterT>::get_nlinsys()
704 : {
705 : return (int)info[9];
706 : }
707 :
708 : template <class fitterT>
709 : std::string levmarInterface<fitterT>::get_reason_string()
710 : {
711 : switch( get_reason_code() )
712 : {
713 : case 1:
714 : return "stopped by small gradient J^T e";
715 : break;
716 : case 2:
717 : return "stopped by small Dp";
718 : break;
719 : case 3:
720 : return "stopped by itmax";
721 : break;
722 : case 4:
723 : return "singular matrix. Restart from current p with increased mu";
724 : break;
725 : case 5:
726 : return "no further error reduction is possible. Restart with increased mu";
727 : break;
728 : case 6:
729 : return "stopped by small ||e||_2";
730 : break;
731 : case 7:
732 : return "stopped by invalid (i.e. NaN or Inf) \"func\" values. This is a user error";
733 : break;
734 : default:
735 : return "unknown reason code";
736 : }
737 : }
738 :
739 : template <class fitterT>
740 : template <typename iosT, char comment>
741 : iosT &levmarInterface<fitterT>::dumpParameters( iosT &ios )
742 : {
743 : // This causes the stream to not output a '\0'
744 : char c[] = { comment, '\0' };
745 :
746 : ios << c << "Current parameters (initial):\n";
747 : for( int i = 0; i < m; i++ )
748 : {
749 : ios << c;
750 : ios << "p[" << i << "] = " << p[i] << " (" << init_p[i] << ")\n";
751 : }
752 :
753 : return ios;
754 : }
755 :
756 : template <class fitterT>
757 : template <char comment>
758 : std::ostream &levmarInterface<fitterT>::dumpParameters()
759 : {
760 : return dumpParameters<std::ostream, comment>( std::cout );
761 : }
762 :
763 : template <class fitterT>
764 : template <typename iosT, char comment>
765 : iosT &levmarInterface<fitterT>::dumpReport( iosT &ios, bool dumpParams )
766 : {
767 : char c[] = { comment, '\0' }; // So a '\0' won't be written to stream.
768 :
769 : ios << c << "--------------------------------------\n";
770 : ios << c << "mx::math::fit::levmarInterface Results \n";
771 : ios << c << "--------------------------------------\n";
772 : if( dumpParams )
773 : dumpParameters<iosT, comment>( ios );
774 : ios << c << "Reason for termination: " << get_reason_string() << "\n";
775 : ios << c << "Initial norm: " << get_initial_norm() << "\n";
776 : ios << c << "Final norm: " << get_final_norm() << "\n";
777 : ios << c << "Number of iterations: " << get_iterations() << "\n";
778 : ios << c << "Function evals: " << get_fevals() << "\n";
779 : ios << c << "Jacobian evals: " << get_jevals() << "\n";
780 : ios << c << "Elapsed time: " << get_deltaT() << " secs\n";
781 : dumpGitStatus<iosT, comment>( ios );
782 : return ios;
783 : }
784 :
785 : template <class fitterT>
786 : template <char comment>
787 : std::ostream &levmarInterface<fitterT>::dumpReport( bool dumpParams )
788 : {
789 : return dumpReport<std::ostream, comment>( std::cout );
790 : }
791 :
792 : /// Test whether a function type has a Jacobian function by testing whether it has a typedef of "hasJacobian"
793 : /** Used for compile-time determination of whether the fitter has a Jacobian.
794 : *
795 : * \ingroup fitting
796 : */
797 : template <typename T>
798 : struct hasJacobian // This was taken directly from the example at
799 : // http://en.wikipedia.org/wiki/Substitution_failure_is_not_an_error
800 : {
801 : // Types "yes" and "no" are guaranteed to have different sizes,
802 : // specifically sizeof(yes) == 1 and sizeof(no) == 2.
803 : typedef char yes[1];
804 : typedef char no[2];
805 :
806 : template <typename fitterT>
807 : static yes &test( typename fitterT::hasJacobian * );
808 :
809 : template <typename>
810 : static no &test( ... );
811 :
812 : /// If hasJacobian<fitterT>::value == true, then fitterT has a Jacobian and the appropriate levmar routines are
813 : /// used. If ::value == false, then the numerical derivatives are calculated.
814 : // If the "sizeof" of the result of calling test<T>(0) would be equal to sizeof(yes),
815 : // the first overload worked and T has a nested type named "hasJacobian".
816 : static const bool value = sizeof( test<T>( 0 ) ) == sizeof( yes );
817 : };
818 :
819 : /// Wrapper for a native array to pass to \ref levmarInterface
820 : /**
821 : * \ingroup fitting
822 : */
823 : template <typename realT>
824 : struct array2Fit
825 : {
826 : realT *data{ 0 }; /// Pointer to the array
827 : size_t nx{ 0 }; /// X dimension of the array
828 : size_t ny{ 0 }; /// Y dimension of the array
829 : };
830 :
831 : } // namespace fit
832 : } // namespace math
833 : } // namespace mx
834 :
835 : #endif // levmarInterface_hpp
|
