<html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>HippoDraw Class Library</title>
<link href="doxygen.css" rel="stylesheet" type="text/css">
<link href="tabs.css" rel="stylesheet" type="text/css">
<link rel="shortcut icon" href="hippoApp.ico">
</head><body>
<!-- Generated by Doxygen 1.6.2 -->
<div class="navigation" id="top">
<div class="tabs">
<ul>
<li><a href="index.html"><span>Main Page</span></a></li>
<li><a href="pages.html"><span>Related Pages</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="annotated.html"><span>Classes</span></a></li>
<li class="current"><a href="files.html"><span>Files</span></a></li>
<li><a href="dirs.html"><span>Directories</span></a></li>
</ul>
</div>
<div class="tabs">
<ul>
<li><a href="files.html"><span>File List</span></a></li>
<li><a href="globals.html"><span>File Members</span></a></li>
</ul>
</div>
<div class="navpath"><a class="el" href="dir_6633238f5107fd0d748d9093e08c459d.html">minimizers</a>
</div>
</div>
<div class="contents">
<h1>BFGSFitter.cxx</h1><a href="BFGSFitter_8cxx.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001
<a name="l00012"></a>00012 <span class="preprocessor">#ifdef _MSC_VER</span>
<a name="l00013"></a>00013 <span class="preprocessor"></span><span class="preprocessor">#include "msdevstudio/MSconfig.h"</span>
<a name="l00014"></a>00014 <span class="preprocessor">#endif</span>
<a name="l00015"></a>00015 <span class="preprocessor"></span>
<a name="l00016"></a>00016 <span class="preprocessor">#ifdef _MSC_VER</span>
<a name="l00017"></a>00017 <span class="preprocessor"></span><span class="preprocessor">#define isnan _isnan</span>
<a name="l00018"></a>00018 <span class="preprocessor"></span><span class="preprocessor">#endif</span>
<a name="l00019"></a>00019 <span class="preprocessor"></span>
<a name="l00020"></a>00020 <span class="comment">//To have isnan.</span>
<a name="l00021"></a>00021 <span class="preprocessor">#ifdef __APPLE__</span>
<a name="l00022"></a>00022 <span class="preprocessor"></span><span class="preprocessor">#include <cstdlib></span>
<a name="l00023"></a>00023 <span class="preprocessor">#define _GLIBCPP_USE_C99 1</span>
<a name="l00024"></a>00024 <span class="preprocessor"></span><span class="preprocessor">#endif</span>
<a name="l00025"></a>00025 <span class="preprocessor"></span>
<a name="l00026"></a>00026 <span class="preprocessor">#include "<a class="code" href="BFGSFitter_8h.html" title="BFGSFitter class interface.">BFGSFitter.h</a>"</span>
<a name="l00027"></a>00027
<a name="l00028"></a>00028 <span class="preprocessor">#include "<a class="code" href="NumLinAlg_8h.html" title="Collection of linear algebra functions.">NumLinAlg.h</a>"</span>
<a name="l00029"></a>00029 <span class="preprocessor">#include "<a class="code" href="StatedFCN_8h.html" title="hippodraw::StatedFCN class interface">StatedFCN.h</a>"</span>
<a name="l00030"></a>00030
<a name="l00031"></a>00031 <span class="comment">// see todo</span>
<a name="l00032"></a>00032 <span class="preprocessor">#include <iostream></span>
<a name="l00033"></a>00033 <span class="keyword">using</span> std::cout;
<a name="l00034"></a>00034 <span class="keyword">using</span> std::endl;
<a name="l00035"></a>00035
<a name="l00036"></a>00036 <span class="preprocessor">#include <cfloat></span>
<a name="l00037"></a>00037 <span class="preprocessor">#include <cstdlib></span>
<a name="l00038"></a>00038 <span class="preprocessor">#include <cassert></span>
<a name="l00039"></a>00039 <span class="preprocessor">#include <cmath></span>
<a name="l00040"></a>00040 <span class="preprocessor">#include <ctime></span>
<a name="l00041"></a>00041
<a name="l00042"></a>00042 <span class="preprocessor">#ifdef __APPLE__</span>
<a name="l00043"></a>00043 <span class="preprocessor"></span><span class="keyword">using</span> std::isnan;
<a name="l00044"></a>00044 <span class="preprocessor">#endif</span>
<a name="l00045"></a>00045 <span class="preprocessor"></span>
<a name="l00046"></a>00046
<a name="l00047"></a>00047 <span class="keyword">using</span> std::pow;
<a name="l00048"></a>00048 <span class="keyword">using</span> std::swap;
<a name="l00049"></a>00049 <span class="keyword">using</span> std::min;
<a name="l00050"></a>00050 <span class="keyword">using</span> std::abs;
<a name="l00051"></a>00051 <span class="keyword">using</span> std::vector;
<a name="l00052"></a>00052 <span class="keyword">using</span> std::vector;
<a name="l00053"></a>00053 <span class="keyword">using</span> std::string;
<a name="l00054"></a>00054 <span class="keyword">using</span> std::map;
<a name="l00055"></a>00055
<a name="l00056"></a>00056 <span class="keyword">using namespace </span>hippodraw;
<a name="l00057"></a>00057
<a name="l00058"></a>00058 <span class="keyword">using namespace </span>Numeric;
<a name="l00059"></a>00059
<a name="l00060"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a0555c889a72414cf283fc0da2ec784bc">00060</a> BFGSFitter::BFGSFitter( <span class="keyword">const</span> <span class="keywordtype">char</span> * name )
<a name="l00061"></a>00061 : <a class="code" href="classhippodraw_1_1Fitter.html" title="The base class for fitters.">Fitter</a> ( name ),
<a name="l00062"></a>00062 m_xinit( 1 ),
<a name="l00063"></a>00063 m_grad_cutoff( 1e-6 ),
<a name="l00064"></a>00064 m_step_cutoff( 1e-6 ),
<a name="l00065"></a>00065 m_proj_cutoff( 1e-6 ),
<a name="l00066"></a>00066 m_c1( 1e-4 ),
<a name="l00067"></a>00067 m_c2( 0.9 ),
<a name="l00068"></a>00068 m_alpha_max( 4 ),
<a name="l00069"></a>00069 m_alpha1( 1 )
<a name="l00070"></a>00070 {
<a name="l00071"></a>00071 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[ <span class="stringliteral">"grad_cutoff"</span> ] = & <a class="code" href="classhippodraw_1_1BFGSFitter.html#a8387240b9424a1270ce2b8fa99d946ec" title="The gradient cut-off parameter.">m_grad_cutoff</a>;
<a name="l00072"></a>00072 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[ <span class="stringliteral">"step_cutoff"</span> ] = & <a class="code" href="classhippodraw_1_1BFGSFitter.html#a0626d59d9cbe05fc96be9119c8b59093" title="The step cut-off parameter.">m_step_cutoff</a>;
<a name="l00073"></a>00073 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[ <span class="stringliteral">"proj_cutoff"</span> ] = & <a class="code" href="classhippodraw_1_1BFGSFitter.html#a21a83b98c0979dc660fe518e2d7c7101" title="The projection cut-off parameter.">m_proj_cutoff</a>;
<a name="l00074"></a>00074 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[ <span class="stringliteral">"c1"</span> ] = & <a class="code" href="classhippodraw_1_1BFGSFitter.html#ad553b78e944ffd23edc383e7eb132db8" title="c1,c2 - constants such that 0 &lt; c1 &lt; c2 &lt; 1 and they ensure that strong...">m_c1</a>;
<a name="l00075"></a>00075 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[ <span class="stringliteral">"c2"</span> ] = & <a class="code" href="classhippodraw_1_1BFGSFitter.html#a916fef9bae3e5d201ee8ba6afc9bdb33">m_c2</a>;
<a name="l00076"></a>00076 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[ <span class="stringliteral">"alpha_max"</span> ] = & <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac6c9defc2c2f41b5f661fb2a80303093" title="Maximum step length to try, suggested value by Nocedal and Wright is alpha_max =...">m_alpha_max</a>;
<a name="l00077"></a>00077 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[ <span class="stringliteral">"alpha1"</span> ] = & <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1bc2b38633aa98e1a99fabdc89636382" title="First step length to try and this must be less than Alpha_max.">m_alpha1</a>;
<a name="l00078"></a>00078 }
<a name="l00079"></a>00079
<a name="l00080"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a53eb174f82d157030b47e1446f597dcd">00080</a> <a class="code" href="classhippodraw_1_1Fitter.html" title="The base class for fitters.">Fitter</a> * <a class="code" href="classhippodraw_1_1BFGSFitter.html#a53eb174f82d157030b47e1446f597dcd" title="Makes a copy of the receiving object.">BFGSFitter::clone</a> ( )<span class="keyword"> const</span>
<a name="l00081"></a>00081 <span class="keyword"></span>{
<a name="l00082"></a>00082 <span class="keywordflow">return</span> <span class="keyword">new</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a0555c889a72414cf283fc0da2ec784bc" title="The constructor taking name of fitter as argument.">BFGSFitter</a> ( *<span class="keyword">this</span> ); <span class="comment">// uses compiler generated copy constructor</span>
<a name="l00083"></a>00083 }
<a name="l00084"></a>00084
<a name="l00087"></a>00087 <span class="keywordtype">bool</span>
<a name="l00088"></a>00088 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a132275ecd9f2fd7e45b34105080ba7bb" title="Main driver routine for BFGS algorithm which has been used in computing the bets...">BFGSFitter::</a>
<a name="l00089"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a132275ecd9f2fd7e45b34105080ba7bb">00089</a> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a132275ecd9f2fd7e45b34105080ba7bb" title="Main driver routine for BFGS algorithm which has been used in computing the bets...">calcBestFit</a>()
<a name="l00090"></a>00090 {
<a name="l00091"></a>00091 <span class="keywordtype">double</span> Alpha_star;
<a name="l00092"></a>00092
<a name="l00093"></a>00093 <span class="comment">// Initialization</span>
<a name="l00094"></a>00094 vector < double > init_params;
<a name="l00095"></a>00095 <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> <a class="code" href="classhippodraw_1_1Fitter.html#a5ff9860d858cb80f4db7e42c08b9df93" title="Fills the vector with the free parameters values.">fillFreeParameters</a> ( init_params );
<a name="l00096"></a>00096 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a7e3fe361ceb3b39af40d35977de4a578" title="Sets the initial value of the iterate, assuming it is given as a vector.">setInitIter</a>( init_params );
<a name="l00097"></a>00097
<a name="l00098"></a>00098 <span class="comment">// Allocate space for various matrices / vectors needed</span>
<a name="l00099"></a>00099 vector< double > xold = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>;
<a name="l00100"></a>00100 vector< double > xnew = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>;
<a name="l00101"></a>00101 vector< double > gk = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a5fadecb4119b5eaa2ef6a36927ba1068" title="The gradient of the objective function.">gradient</a>( xold );
<a name="l00102"></a>00102 vector< double > gkPlusUn = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a5fadecb4119b5eaa2ef6a36927ba1068" title="The gradient of the objective function.">gradient</a>( xnew );
<a name="l00103"></a>00103
<a name="l00104"></a>00104 vector< double > pk( <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size() );
<a name="l00105"></a>00105 vector< double > s( <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size() );
<a name="l00106"></a>00106 vector< double > y( <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size() );
<a name="l00107"></a>00107
<a name="l00108"></a>00108 <span class="comment">// The standard quasi newton initialization </span>
<a name="l00109"></a>00109 <a class="code" href="namespacehippodraw_1_1Numeric.html#a3c1c282c14b12b9f982f4f3f359fde67" title="Creates an n x n identity matrix for M.">eye</a> ( <a class="code" href="classhippodraw_1_1BFGSFitter.html#a68f9abfe3d088353c41bf490c1ef5339" title="The inverse of the quasi-Hessian.">m_M</a>, <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size() );
<a name="l00110"></a>00110 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a68f9abfe3d088353c41bf490c1ef5339" title="The inverse of the quasi-Hessian.">m_M</a> = ( 1.0 / <a class="code" href="namespacehippodraw_1_1Numeric.html#a2fc63b8199da4a6f45b4a9590fbab38f" title="Computes the two norm of the vector.">norm</a>( gk ) ) * m_M;
<a name="l00111"></a>00111
<a name="l00112"></a>00112 vector< vector< double > > t1 , t2;
<a name="l00113"></a>00113 <a class="code" href="namespacehippodraw_1_1Numeric.html#a3c1c282c14b12b9f982f4f3f359fde67" title="Creates an n x n identity matrix for M.">eye</a>( t1, <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size() );
<a name="l00114"></a>00114 <a class="code" href="namespacehippodraw_1_1Numeric.html#a3c1c282c14b12b9f982f4f3f359fde67" title="Creates an n x n identity matrix for M.">eye</a>( t2, <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size() );
<a name="l00115"></a>00115
<a name="l00116"></a>00116 <span class="keywordtype">double</span> fx = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( xnew );
<a name="l00117"></a>00117 <span class="keywordtype">double</span> fxold = fx;
<a name="l00118"></a>00118 <span class="keywordflow">for</span>( <span class="keywordtype">int</span> k = 1; k <= <a class="code" href="classhippodraw_1_1Fitter.html#a35b1e829c1c13a7aef0c1b4cc1895c34" title="The maximum number of iterations allowed in attempting the fit.">m_max_iterations</a>; k++ )
<a name="l00119"></a>00119 {
<a name="l00120"></a>00120 <span class="comment">// Update the iterate.</span>
<a name="l00121"></a>00121 gk = gkPlusUn;
<a name="l00122"></a>00122 pk = m_M * (-gk);
<a name="l00123"></a>00123 Alpha_star = <a class="code" href="classhippodraw_1_1BFGSFitter.html#ae62893d4bdeceb6eaf764a4ef31458fb" title="Computes a step satisfying the Wolfe conditions.">wolfeStep</a>( xold, pk );
<a name="l00124"></a>00124
<a name="l00125"></a>00125 <span class="keywordflow">do</span>
<a name="l00126"></a>00126 {
<a name="l00127"></a>00127 xnew = xold + Alpha_star * pk;
<a name="l00128"></a>00128 fx = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( xnew );
<a name="l00129"></a>00129 Alpha_star /= 3.0;
<a name="l00130"></a>00130 }
<a name="l00131"></a>00131 <span class="keywordflow">while</span>( isnan( fx ) );
<a name="l00132"></a>00132
<a name="l00133"></a>00133 gkPlusUn = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a5fadecb4119b5eaa2ef6a36927ba1068" title="The gradient of the objective function.">gradient</a>( xnew );
<a name="l00134"></a>00134
<a name="l00135"></a>00135 s = xnew - xold;
<a name="l00136"></a>00136 y = gkPlusUn - gk;
<a name="l00137"></a>00137
<a name="l00138"></a>00138 <span class="keywordtype">double</span> ys = <a class="code" href="namespacehippodraw_1_1Numeric.html#a329ce3200d86ee77e2b9721db28388ee" title="Computes the dot or the inner product of two vectors(i.e.">innerProduct</a>( y, s );
<a name="l00139"></a>00139 <span class="keywordtype">double</span> yy = <a class="code" href="namespacehippodraw_1_1Numeric.html#a2fc63b8199da4a6f45b4a9590fbab38f" title="Computes the two norm of the vector.">norm</a>( y );
<a name="l00140"></a>00140 <span class="keywordtype">double</span> ss = <a class="code" href="namespacehippodraw_1_1Numeric.html#a2fc63b8199da4a6f45b4a9590fbab38f" title="Computes the two norm of the vector.">norm</a>( s );
<a name="l00141"></a>00141
<a name="l00142"></a>00142 <span class="comment">// Termination criteria</span>
<a name="l00143"></a>00143 <span class="keywordflow">if</span>( ( abs( ys ) < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a21a83b98c0979dc660fe518e2d7c7101" title="The projection cut-off parameter.">m_proj_cutoff</a> ) ||
<a name="l00144"></a>00144 ( abs( Alpha_star ) < DBL_EPSILON ) ||
<a name="l00145"></a>00145 ( ss < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a0626d59d9cbe05fc96be9119c8b59093" title="The step cut-off parameter.">m_step_cutoff</a> ) ||
<a name="l00146"></a>00146 ( yy < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a8387240b9424a1270ce2b8fa99d946ec" title="The gradient cut-off parameter.">m_grad_cutoff</a> ) ||
<a name="l00147"></a>00147 ( fx >= fxold ) )
<a name="l00148"></a>00148 <span class="keywordflow">break</span>;
<a name="l00149"></a>00149
<a name="l00150"></a>00150 <span class="comment">// DFP Update of inverse of approximate hessian.</span>
<a name="l00151"></a>00151 <span class="comment">// M = M-(s*y'*M+M*y*s')/(y'*s)+(1+(y'*M*y)/(y'*s))*(s*s')/(y'*s);</span>
<a name="l00152"></a>00152 <span class="keywordtype">double</span> temp = ( 1 + <a class="code" href="namespacehippodraw_1_1Numeric.html#a329ce3200d86ee77e2b9721db28388ee" title="Computes the dot or the inner product of two vectors(i.e.">innerProduct</a>( y, m_M * y ) / ys ) / ys;
<a name="l00153"></a>00153
<a name="l00154"></a>00154 t1 = ( <a class="code" href="namespacehippodraw_1_1Numeric.html#af9cb017e1c9362831eccadf6bceac017" title="Computes the outer product of two vectors (i.e.">outerProduct</a>(s, y) * m_M)/ys + ( m_M * <a class="code" href="namespacehippodraw_1_1Numeric.html#af9cb017e1c9362831eccadf6bceac017" title="Computes the outer product of two vectors (i.e.">outerProduct</a>(y , s) ) / ys;
<a name="l00155"></a>00155 t2 = temp * <a class="code" href="namespacehippodraw_1_1Numeric.html#af9cb017e1c9362831eccadf6bceac017" title="Computes the outer product of two vectors (i.e.">outerProduct</a>( s, s );
<a name="l00156"></a>00156 m_M = m_M - t1 + t2;
<a name="l00157"></a>00157
<a name="l00158"></a>00158 <span class="comment">// one pass of the loop is over so refresh</span>
<a name="l00159"></a>00159 xold = xnew;
<a name="l00160"></a>00160 fxold = fx;
<a name="l00161"></a>00161 }
<a name="l00162"></a>00162
<a name="l00163"></a>00163 <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> setFreeParameters ( xold );
<a name="l00164"></a>00164 <span class="comment">//write( xold );</span>
<a name="l00165"></a>00165
<a name="l00166"></a>00166 <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00167"></a>00167 }
<a name="l00168"></a>00168
<a name="l00171"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#ae62893d4bdeceb6eaf764a4ef31458fb">00171</a> <span class="keywordtype">double</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#ae62893d4bdeceb6eaf764a4ef31458fb" title="Computes a step satisfying the Wolfe conditions.">BFGSFitter::wolfeStep</a>( <span class="keyword">const</span> std::vector< double >& x0,
<a name="l00172"></a>00172 <span class="keyword">const</span> std::vector< double >& p )<span class="keyword"> const</span>
<a name="l00173"></a>00173 <span class="keyword"></span>{
<a name="l00174"></a>00174 <span class="keywordtype">double</span> step_fac = sqrt(2.0); <span class="comment">// Geometric step factor; must be > 1</span>
<a name="l00175"></a>00175
<a name="l00176"></a>00176 <span class="keywordtype">double</span> phi0 = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( x0 ); <span class="comment">// Function value at the initial point</span>
<a name="l00177"></a>00177 <span class="keywordtype">double</span> dphi0 = <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac70817d96bec5f159f849a35eca80317" title="Efficient computation of gradient of the objective function with a vector p.">gradp</a>( x0, p ); <span class="comment">// innner product of gradient and p </span>
<a name="l00178"></a>00178
<a name="l00179"></a>00179 <span class="comment">// dphi0 >= 0 means you are ascending. To avoid it set step = 0 and</span>
<a name="l00180"></a>00180 <span class="comment">// terminate the iterations ( maybe after giving out a warning )</span>
<a name="l00181"></a>00181 <span class="comment">// Algorithm ensures that such case will not arise but curse of finite</span>
<a name="l00182"></a>00182 <span class="comment">// precision mathematics lead to some very small positive dphi0.</span>
<a name="l00183"></a>00183 <span class="keywordflow">if</span> ( dphi0 >= 0 )
<a name="l00184"></a>00184 <span class="keywordflow">return</span> 0.0;
<a name="l00185"></a>00185
<a name="l00186"></a>00186 <span class="keywordtype">double</span> Alpha0 = 0;
<a name="l00187"></a>00187 <span class="keywordtype">double</span> Alphai = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1bc2b38633aa98e1a99fabdc89636382" title="First step length to try and this must be less than Alpha_max.">m_alpha1</a>;
<a name="l00188"></a>00188 <span class="keywordtype">double</span> Alphaim = Alpha0;
<a name="l00189"></a>00189 <span class="keywordtype">double</span> phiim = phi0;
<a name="l00190"></a>00190 <span class="keywordtype">int</span> i = 1;
<a name="l00191"></a>00191 <span class="keywordtype">int</span> done = 0;
<a name="l00192"></a>00192 <span class="keywordtype">int</span> cnt = 0;
<a name="l00193"></a>00193
<a name="l00194"></a>00194 <span class="keywordtype">double</span> phii, dphii;
<a name="l00195"></a>00195
<a name="l00196"></a>00196 <span class="keywordflow">while</span> ( !done )
<a name="l00197"></a>00197 {
<a name="l00198"></a>00198 phii = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( x0 + Alphai * p );
<a name="l00199"></a>00199
<a name="l00200"></a>00200 <span class="keywordflow">if</span> ( (phii > (phi0 + <a class="code" href="classhippodraw_1_1BFGSFitter.html#ad553b78e944ffd23edc383e7eb132db8" title="c1,c2 - constants such that 0 &lt; c1 &lt; c2 &lt; 1 and they ensure that strong...">m_c1</a> * Alphai * dphi0)) ||
<a name="l00201"></a>00201 ( (phii >= phiim) && ( i > 1)) )
<a name="l00202"></a>00202 <span class="keywordflow">return</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#adfacd0c34f0889cef880bf7f82394751" title="A function which helps out Wolfe in deciding the step length.">zoom</a>( x0, p, phi0, dphi0, Alphaim, Alphai );
<a name="l00203"></a>00203
<a name="l00204"></a>00204 dphii = <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac70817d96bec5f159f849a35eca80317" title="Efficient computation of gradient of the objective function with a vector p.">gradp</a>( x0 + Alphai * p , p );
<a name="l00205"></a>00205
<a name="l00206"></a>00206 <span class="keywordflow">if</span> ( abs( dphii ) <= -<a class="code" href="classhippodraw_1_1BFGSFitter.html#a916fef9bae3e5d201ee8ba6afc9bdb33">m_c2</a> * dphi0 )
<a name="l00207"></a>00207 <span class="keywordflow">return</span> Alphai;
<a name="l00208"></a>00208
<a name="l00209"></a>00209 <span class="keywordflow">if</span> (dphii >= 0)
<a name="l00210"></a>00210 <span class="keywordflow">return</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#adfacd0c34f0889cef880bf7f82394751" title="A function which helps out Wolfe in deciding the step length.">zoom</a>( x0, p, phi0, dphi0, Alphai, Alphaim );
<a name="l00211"></a>00211
<a name="l00212"></a>00212 <span class="comment">// Choose new Alphai in (Alphai, Alpha_max)</span>
<a name="l00213"></a>00213 Alphaim = Alphai;
<a name="l00214"></a>00214 phiim = phii;
<a name="l00215"></a>00215 Alphai = min( step_fac * Alphai, <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac6c9defc2c2f41b5f661fb2a80303093" title="Maximum step length to try, suggested value by Nocedal and Wright is alpha_max =...">m_alpha_max</a> );
<a name="l00216"></a>00216
<a name="l00217"></a>00217 <span class="keywordflow">if</span> ( Alphai == <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac6c9defc2c2f41b5f661fb2a80303093" title="Maximum step length to try, suggested value by Nocedal and Wright is alpha_max =...">m_alpha_max</a> )
<a name="l00218"></a>00218 cnt += 1;
<a name="l00219"></a>00219
<a name="l00220"></a>00220 <span class="keywordflow">if</span> (cnt > 1)
<a name="l00221"></a>00221 {
<a name="l00222"></a>00222 <span class="comment">// cout << "WARNING: Unable to bracket a strong Wolfe pt. in [ "</span>
<a name="l00223"></a>00223 <span class="comment">// << m_alpha1 << ", " << m_alpha_max << " ]" << endl;</span>
<a name="l00224"></a>00224 <span class="keywordflow">return</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac6c9defc2c2f41b5f661fb2a80303093" title="Maximum step length to try, suggested value by Nocedal and Wright is alpha_max =...">m_alpha_max</a>;
<a name="l00225"></a>00225 }
<a name="l00226"></a>00226
<a name="l00227"></a>00227 i += 1;
<a name="l00228"></a>00228 <span class="comment">//cout << " Wolfe Iteration: " << i << endl;</span>
<a name="l00229"></a>00229 }
<a name="l00230"></a>00230
<a name="l00231"></a>00231 <span class="keywordflow">return</span> 0.0;
<a name="l00232"></a>00232 }
<a name="l00233"></a>00233
<a name="l00234"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#adfacd0c34f0889cef880bf7f82394751">00234</a> <span class="keywordtype">double</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#adfacd0c34f0889cef880bf7f82394751" title="A function which helps out Wolfe in deciding the step length.">BFGSFitter::zoom</a>( <span class="keyword">const</span> std::vector< double >& x0,
<a name="l00235"></a>00235 <span class="keyword">const</span> std::vector< double >& p,
<a name="l00236"></a>00236 <span class="keywordtype">double</span> phi0, <span class="keywordtype">double</span> dphi0,
<a name="l00237"></a>00237 <span class="keywordtype">double</span> Alphalo, <span class="keywordtype">double</span> Alphahi )<span class="keyword"> const</span>
<a name="l00238"></a>00238 <span class="keyword"></span>{
<a name="l00239"></a>00239 <span class="keywordtype">int</span> MaxIter = 20;
<a name="l00240"></a>00240 <span class="keywordtype">int</span> iter = 0;
<a name="l00241"></a>00241 <span class="keywordtype">int</span> done = 0;
<a name="l00242"></a>00242
<a name="l00243"></a>00243 <span class="keywordtype">double</span> philo, phij;
<a name="l00244"></a>00244 <span class="keywordtype">double</span> dphij;
<a name="l00245"></a>00245 <span class="keywordtype">double</span> Alphaj;
<a name="l00246"></a>00246 <span class="keywordtype">double</span> Alpha_star = 0.0;
<a name="l00247"></a>00247
<a name="l00248"></a>00248 <span class="keywordflow">while</span> ( !done && iter < MaxIter )
<a name="l00249"></a>00249 {
<a name="l00250"></a>00250 iter += 1;
<a name="l00251"></a>00251
<a name="l00252"></a>00252 philo = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( x0 + Alphalo * p );
<a name="l00253"></a>00253
<a name="l00254"></a>00254 <span class="comment">// Find a trial step length Alphaj between Alphalo and Alphahi</span>
<a name="l00255"></a>00255 Alphaj = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a31514e4ebc5aaf72ff8f44a35a122120" title="A cubic interpolation routine.">interpolate</a>( x0, p, Alphalo, Alphahi );
<a name="l00256"></a>00256
<a name="l00257"></a>00257 <span class="comment">// Evalaute phi( Alphaj )</span>
<a name="l00258"></a>00258 phij = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( x0 + Alphaj * p );
<a name="l00259"></a>00259
<a name="l00260"></a>00260 <span class="keywordflow">if</span>( (phij > phi0 + <a class="code" href="classhippodraw_1_1BFGSFitter.html#ad553b78e944ffd23edc383e7eb132db8" title="c1,c2 - constants such that 0 &lt; c1 &lt; c2 &lt; 1 and they ensure that strong...">m_c1</a> * Alphaj * dphi0) || (phij >= philo) )
<a name="l00261"></a>00261 Alphahi = Alphaj;
<a name="l00262"></a>00262 <span class="keywordflow">else</span>
<a name="l00263"></a>00263 {
<a name="l00264"></a>00264 dphij = <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac70817d96bec5f159f849a35eca80317" title="Efficient computation of gradient of the objective function with a vector p.">gradp</a>( x0 + Alphaj * p , p );
<a name="l00265"></a>00265
<a name="l00266"></a>00266 <span class="keywordflow">if</span> ( abs( dphij ) <= -<a class="code" href="classhippodraw_1_1BFGSFitter.html#a916fef9bae3e5d201ee8ba6afc9bdb33">m_c2</a> * dphi0 )
<a name="l00267"></a>00267 Alpha_star = Alphaj;
<a name="l00268"></a>00268 <span class="keywordflow">return</span> Alpha_star;
<a name="l00269"></a>00269
<a name="l00270"></a>00270 <span class="keywordflow">if</span> ( dphij * ( Alphahi - Alphalo ) >= 0 )
<a name="l00271"></a>00271 Alphahi = Alphalo;
<a name="l00272"></a>00272
<a name="l00273"></a>00273 Alphalo = Alphaj;
<a name="l00274"></a>00274 }
<a name="l00275"></a>00275
<a name="l00276"></a>00276 }
<a name="l00277"></a>00277
<a name="l00278"></a>00278 <span class="comment">// If above loop fails take the mid-point</span>
<a name="l00279"></a>00279 <span class="keywordflow">if</span> (iter == MaxIter)
<a name="l00280"></a>00280 Alpha_star = 0.5 * ( Alphahi + Alphalo );
<a name="l00281"></a>00281
<a name="l00282"></a>00282 <span class="keywordflow">return</span> Alpha_star;
<a name="l00283"></a>00283
<a name="l00284"></a>00284 }
<a name="l00285"></a>00285
<a name="l00286"></a>00286
<a name="l00287"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a31514e4ebc5aaf72ff8f44a35a122120">00287</a> <span class="keywordtype">double</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a31514e4ebc5aaf72ff8f44a35a122120" title="A cubic interpolation routine.">BFGSFitter::interpolate</a>( <span class="keyword">const</span> std::vector< double >& x0,
<a name="l00288"></a>00288 <span class="keyword">const</span> std::vector< double >& p,
<a name="l00289"></a>00289 <span class="keywordtype">double</span> Alphaim,
<a name="l00290"></a>00290 <span class="keywordtype">double</span> Alphai )<span class="keyword"> const</span>
<a name="l00291"></a>00291 <span class="keyword"></span>{
<a name="l00292"></a>00292
<a name="l00293"></a>00293 <span class="keywordflow">if</span> ( Alphaim > Alphai )
<a name="l00294"></a>00294 swap( Alphaim, Alphai);
<a name="l00295"></a>00295
<a name="l00296"></a>00296 <span class="keywordtype">double</span> phiim = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( x0 + Alphaim * p );
<a name="l00297"></a>00297 <span class="keywordtype">double</span> phii = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">function</a>( x0 + Alphai * p );
<a name="l00298"></a>00298
<a name="l00299"></a>00299 <span class="keywordtype">double</span> dphiim = <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac70817d96bec5f159f849a35eca80317" title="Efficient computation of gradient of the objective function with a vector p.">gradp</a>( x0 + Alphaim * p, p );
<a name="l00300"></a>00300 <span class="keywordtype">double</span> dphii = <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac70817d96bec5f159f849a35eca80317" title="Efficient computation of gradient of the objective function with a vector p.">gradp</a>( x0 + Alphai * p, p );
<a name="l00301"></a>00301
<a name="l00302"></a>00302 <span class="keywordtype">double</span> d1 = dphiim + dphii - 3 * ( (phiim - phii)/(Alphaim - Alphai) );
<a name="l00303"></a>00303 <span class="keywordtype">double</span> d2 = sqrt( d1 * d1 - dphiim * dphii);
<a name="l00304"></a>00304
<a name="l00305"></a>00305 <span class="keywordtype">double</span> Alphaip = Alphai - (Alphai - Alphaim) *
<a name="l00306"></a>00306 ( (dphiim + d2 - d1) / (dphii - dphiim + 2 * d2) );
<a name="l00307"></a>00307
<a name="l00308"></a>00308 <span class="keywordtype">double</span> lth = abs(Alphai - Alphaim);
<a name="l00309"></a>00309
<a name="l00310"></a>00310 <span class="keywordflow">if</span>( abs(Alphaip - Alphai) < 0.05 * lth ||
<a name="l00311"></a>00311 abs(Alphaip - Alphaim) < 0.05 * lth ||
<a name="l00312"></a>00312 Alphaip < Alphaim ||
<a name="l00313"></a>00313 Alphaip > Alphai )
<a name="l00314"></a>00314 Alphaip = 0.5 * (Alphai + Alphaim);
<a name="l00315"></a>00315
<a name="l00316"></a>00316 <span class="keywordflow">return</span> Alphaip;
<a name="l00317"></a>00317 }
<a name="l00318"></a>00318
<a name="l00319"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320">00319</a> <span class="keywordtype">double</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1618f215be470523594cc47d286f1320" title="The objective function.">BFGSFitter::function</a>( <span class="keyword">const</span> std::vector< double > & u )<span class="keyword"> const</span>
<a name="l00320"></a>00320 <span class="keyword"></span>{
<a name="l00321"></a>00321 <span class="comment">// Gets the value of the objective function from pvfcn</span>
<a name="l00322"></a>00322 <span class="comment">// first you have to convert the vector to vector</span>
<a name="l00323"></a>00323 vector< double > x( u.size() );
<a name="l00324"></a>00324
<a name="l00325"></a>00325 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < u.size(); i++ )
<a name="l00326"></a>00326 x[i] = u[i];
<a name="l00327"></a>00327
<a name="l00328"></a>00328 <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> setFreeParameters ( x );
<a name="l00329"></a>00329
<a name="l00330"></a>00330 <span class="keywordflow">return</span> <a class="code" href="classhippodraw_1_1Fitter.html#a35d807e079d04d56ca61b604d49deae8" title="Calculates the value of the objective function at the current set of parameters.">objectiveValue</a>();
<a name="l00331"></a>00331
<a name="l00332"></a>00332 <span class="comment">// Following are a few standard test functions</span>
<a name="l00333"></a>00333 <span class="comment">// which were used to test this minimizer. </span>
<a name="l00334"></a>00334
<a name="l00335"></a>00335 <span class="comment">//double x = u[0]; double y = u[1];</span>
<a name="l00336"></a>00336
<a name="l00337"></a>00337 <span class="comment">// Rosenbrock's function</span>
<a name="l00338"></a>00338 <span class="comment">//return 100 * (y - x * x) * (y - x * x) + (1 - x) * (1 - x);</span>
<a name="l00339"></a>00339
<a name="l00340"></a>00340 <span class="comment">// Freudenstein and Roth's Function</span>
<a name="l00341"></a>00341 <span class="comment">//return pow(-13 + x + ((5 - y)*y - 2 )*y, 2) + </span>
<a name="l00342"></a>00342 <span class="comment">// pow(-29 + x + ((y + 1)*y - 14)*y, 2);</span>
<a name="l00343"></a>00343
<a name="l00344"></a>00344 <span class="comment">// Beale Function </span>
<a name="l00345"></a>00345 <span class="comment">//return pow(1.5 - x*(1 - y ), 2) + </span>
<a name="l00346"></a>00346 <span class="comment">// pow(2.25 - x*(1 - y*y ), 2) + </span>
<a name="l00347"></a>00347 <span class="comment">// pow(2.625 - x*(1 - y*y*y ), 2);</span>
<a name="l00348"></a>00348 }
<a name="l00349"></a>00349
<a name="l00350"></a>00350 vector< double >
<a name="l00351"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a5fadecb4119b5eaa2ef6a36927ba1068">00351</a> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a5fadecb4119b5eaa2ef6a36927ba1068" title="The gradient of the objective function.">BFGSFitter::gradient</a>( <span class="keyword">const</span> std::vector< double > & u )<span class="keyword"> const</span>
<a name="l00352"></a>00352 <span class="keyword"></span>{
<a name="l00353"></a>00353 <span class="keywordtype">double</span> h = 1e-5;
<a name="l00354"></a>00354 vector< double > x( u.size(), 0.0 );
<a name="l00355"></a>00355 vector< double > xph( u.size(), 0.0 );
<a name="l00356"></a>00356
<a name="l00357"></a>00357 vector< double > g( u.size() );
<a name="l00358"></a>00358
<a name="l00359"></a>00359 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < u.size(); i++ )
<a name="l00360"></a>00360 x[i] = u[i];
<a name="l00361"></a>00361
<a name="l00362"></a>00362 <span class="comment">// Calculating the gradient by finite differencing</span>
<a name="l00363"></a>00363 <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> setFreeParameters ( x );
<a name="l00364"></a>00364 <span class="keywordtype">double</span> fx = <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> <a class="code" href="classhippodraw_1_1Fitter.html#a35d807e079d04d56ca61b604d49deae8" title="Calculates the value of the objective function at the current set of parameters.">objectiveValue</a> ();
<a name="l00365"></a>00365 <span class="keywordtype">double</span> fxph = 0.0;
<a name="l00366"></a>00366 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < u.size(); i++ )
<a name="l00367"></a>00367 {
<a name="l00368"></a>00368 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < u.size(); j++ ) {
<a name="l00369"></a>00369 xph[j] = ( i == j ) ? ( x[j] + h ) : x[j];
<a name="l00370"></a>00370 }
<a name="l00371"></a>00371 <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> setFreeParameters ( xph );
<a name="l00372"></a>00372 fxph = <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> <a class="code" href="classhippodraw_1_1Fitter.html#a35d807e079d04d56ca61b604d49deae8" title="Calculates the value of the objective function at the current set of parameters.">objectiveValue</a> ();
<a name="l00373"></a>00373 g[i] = ( fxph - fx ) / h;
<a name="l00374"></a>00374 }
<a name="l00375"></a>00375
<a name="l00376"></a>00376 <span class="comment">// Following are a few gradients of standard test functions</span>
<a name="l00377"></a>00377 <span class="comment">// which were used to test this minimizer.</span>
<a name="l00378"></a>00378
<a name="l00379"></a>00379 <span class="comment">/* double x = u[0]; double y = u[1];</span>
<a name="l00380"></a>00380 <span class="comment"> vector< double > g(2);*/</span>
<a name="l00381"></a>00381
<a name="l00382"></a>00382 <span class="comment">// Gradient of Rosenbrock's function</span>
<a name="l00383"></a>00383 <span class="comment">/* g[0] = -400 * x * ( y - x * x ) - 2 * ( 1 - x );</span>
<a name="l00384"></a>00384 <span class="comment"> g[1] = 200 * ( y - x * x );*/</span>
<a name="l00385"></a>00385
<a name="l00386"></a>00386 <span class="comment">// Gradient of Freudenstein and Roth's Function</span>
<a name="l00387"></a>00387 <span class="comment">//g[0] = 2*(y*(y*(y+1)-14)+ x - 29) +</span>
<a name="l00388"></a>00388 <span class="comment">// 2*(y*((5-y)*y-2) + x - 13);</span>
<a name="l00389"></a>00389 <span class="comment">//g[1] = 2*(y*(2*y+1) + y*(y+1) - 14) * (y*(y*(y+1) - 14) + x - 29) + </span>
<a name="l00390"></a>00390 <span class="comment">// 2*((5-y)*y + (5-2*y)*y - 2) * (y*((5-y)*y - 2) + x - 13);</span>
<a name="l00391"></a>00391
<a name="l00392"></a>00392 <span class="comment">// Gradient of Beale Function</span>
<a name="l00393"></a>00393 <span class="comment">//g[0]= 2*(2.625 - x * (1 - y*y*y) ) * (y*y*y - 1)+</span>
<a name="l00394"></a>00394 <span class="comment">// 2*(2.25 - x * (1 - y*y) ) * (y*y - 1)+</span>
<a name="l00395"></a>00395 <span class="comment">// 2*(1.5 - x * (1 - y ) ) * (y - 1);</span>
<a name="l00396"></a>00396 <span class="comment">//g[1]= 6*x*y*y*(2.625 - x * (1 - y*y*y)) + </span>
<a name="l00397"></a>00397 <span class="comment">// 4*x*y *(2.25 - x * (1 - y*y)) + </span>
<a name="l00398"></a>00398 <span class="comment">// 2*x *(1.5 - x * (1 - y ));</span>
<a name="l00399"></a>00399
<a name="l00400"></a>00400 <span class="keywordflow">return</span> g;
<a name="l00401"></a>00401 }
<a name="l00402"></a>00402
<a name="l00403"></a>00403
<a name="l00404"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#ac70817d96bec5f159f849a35eca80317">00404</a> <span class="keywordtype">double</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#ac70817d96bec5f159f849a35eca80317" title="Efficient computation of gradient of the objective function with a vector p.">BFGSFitter::gradp</a>( <span class="keyword">const</span> std::vector< double > & u,
<a name="l00405"></a>00405 <span class="keyword">const</span> std::vector< double > & p )<span class="keyword"> const</span>
<a name="l00406"></a>00406 <span class="keyword"></span>{
<a name="l00407"></a>00407 <span class="keywordtype">double</span> h = 1e-5;
<a name="l00408"></a>00408 vector< double > x( u.size() );
<a name="l00409"></a>00409
<a name="l00410"></a>00410 <span class="comment">// Calculating the gradient in direction of p by finite differencing</span>
<a name="l00411"></a>00411 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < u.size(); i++ ) {
<a name="l00412"></a>00412 x[i] = u[i];
<a name="l00413"></a>00413 }
<a name="l00414"></a>00414 <span class="comment">// double fx = m_fcn -> operator()( x );</span>
<a name="l00415"></a>00415 <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> setFreeParameters ( x );
<a name="l00416"></a>00416 <span class="keywordtype">double</span> fx = <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> <a class="code" href="classhippodraw_1_1Fitter.html#a35d807e079d04d56ca61b604d49deae8" title="Calculates the value of the objective function at the current set of parameters.">objectiveValue</a> ( );
<a name="l00417"></a>00417
<a name="l00418"></a>00418 <span class="keywordflow">for</span> ( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < u.size(); i++ ) {
<a name="l00419"></a>00419 x[i] += h * p[i] ;
<a name="l00420"></a>00420 }
<a name="l00421"></a>00421 <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> setFreeParameters ( x );
<a name="l00422"></a>00422 <span class="keywordtype">double</span> fxph = <a class="code" href="classhippodraw_1_1Fitter.html#a1f88348dd4b35c18128ec042fea57ee2" title="The objective function.">m_fcn</a> -> <a class="code" href="classhippodraw_1_1Fitter.html#a35d807e079d04d56ca61b604d49deae8" title="Calculates the value of the objective function at the current set of parameters.">objectiveValue</a> ();
<a name="l00423"></a>00423
<a name="l00424"></a>00424 <span class="keywordflow">return</span> ( fxph - fx ) / h;
<a name="l00425"></a>00425
<a name="l00426"></a>00426 <span class="comment">// Following are a directional derivative of standard test functions</span>
<a name="l00427"></a>00427 <span class="comment">// which were used to test this minimizer.</span>
<a name="l00428"></a>00428
<a name="l00429"></a>00429 <span class="comment">/* double x = u[0]; double y = u[1];</span>
<a name="l00430"></a>00430 <span class="comment"> vector< double > g(2);*/</span>
<a name="l00431"></a>00431
<a name="l00432"></a>00432 <span class="comment">// Gradient of Rosenbrock's function</span>
<a name="l00433"></a>00433 <span class="comment">/* g[0] = -400 * x * ( y - x * x ) - 2 * ( 1 - x );</span>
<a name="l00434"></a>00434 <span class="comment"> g[1] = 200 * ( y - x * x );*/</span>
<a name="l00435"></a>00435
<a name="l00436"></a>00436 <span class="comment">// Gradient of Freudenstein and Roth's Function</span>
<a name="l00437"></a>00437 <span class="comment">//g[0] = 2*(y*(y*(y+1)-14)+ x - 29) +</span>
<a name="l00438"></a>00438 <span class="comment">// 2*(y*((5-y)*y-2) + x - 13);</span>
<a name="l00439"></a>00439 <span class="comment">//g[1] = 2*(y*(2*y+1) + y*(y+1) - 14) * (y*(y*(y+1) - 14) + x - 29) + </span>
<a name="l00440"></a>00440 <span class="comment">// 2*((5-y)*y + (5-2*y)*y - 2) * (y*((5-y)*y - 2) + x - 13);</span>
<a name="l00441"></a>00441
<a name="l00442"></a>00442 <span class="comment">// Gradient of Beale Function</span>
<a name="l00443"></a>00443 <span class="comment">//g[0]= 2*(2.625 - x * (1 - y*y*y) ) * (y*y*y - 1)+</span>
<a name="l00444"></a>00444 <span class="comment">// 2*(2.25 - x * (1 - y*y) ) * (y*y - 1)+</span>
<a name="l00445"></a>00445 <span class="comment">// 2*(1.5 - x * (1 - y ) ) * (y - 1);</span>
<a name="l00446"></a>00446 <span class="comment">//g[1]= 6*x*y*y*(2.625 - x * (1 - y*y*y)) + </span>
<a name="l00447"></a>00447 <span class="comment">// 4*x*y *(2.25 - x * (1 - y*y)) + </span>
<a name="l00448"></a>00448 <span class="comment">// 2*x *(1.5 - x * (1 - y ));</span>
<a name="l00449"></a>00449
<a name="l00450"></a>00450 <span class="comment">// return g[0] * p[0] + g[1] * p[1];</span>
<a name="l00451"></a>00451 }
<a name="l00452"></a>00452
<a name="l00453"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#add5d43beb0cf4b646200adb03208d6e9">00453</a> <span class="keyword">const</span> vector< double > & <a class="code" href="classhippodraw_1_1BFGSFitter.html#add5d43beb0cf4b646200adb03208d6e9" title="Returns the initial value of the iterate.">BFGSFitter::initIter</a>()<span class="keyword"> const</span>
<a name="l00454"></a>00454 <span class="keyword"></span>{
<a name="l00455"></a>00455 <span class="keywordflow">return</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>;
<a name="l00456"></a>00456 }
<a name="l00457"></a>00457
<a name="l00458"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a7e3fe361ceb3b39af40d35977de4a578">00458</a> <span class="keywordtype">int</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a7e3fe361ceb3b39af40d35977de4a578" title="Sets the initial value of the iterate, assuming it is given as a vector.">BFGSFitter::setInitIter</a>( <span class="keyword">const</span> std::vector< double > & xinit )
<a name="l00459"></a>00459 {
<a name="l00460"></a>00460 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.resize( xinit.size() );
<a name="l00461"></a>00461
<a name="l00462"></a>00462 <span class="comment">// Provide a random perturbations to the initial value</span>
<a name="l00463"></a>00463 <span class="comment">//srand( (unsigned) time( NULL ) );</span>
<a name="l00464"></a>00464 <span class="comment">//for( unsigned int i = 0; i < xinit.size(); i++ )</span>
<a name="l00465"></a>00465 <span class="comment">//m_xinit[i] = xinit[i] * (1 + 0.025 * ( 0.5 - rand() / ( RAND_MAX + 1.0 )));</span>
<a name="l00466"></a>00466
<a name="l00467"></a>00467 <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a> = xinit;
<a name="l00468"></a>00468
<a name="l00469"></a>00469 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00470"></a>00470 }
<a name="l00471"></a>00471
<a name="l00472"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a3ca477145d995d896e7cdce3d23a8c0d">00472</a> <span class="keywordtype">int</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a3ca477145d995d896e7cdce3d23a8c0d" title="Calculates the covariance matrix.">BFGSFitter::calcCovariance</a> ( std::vector < std::vector < double > >& cov )
<a name="l00473"></a>00473 {
<a name="l00474"></a>00474 cov.resize( <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size() );
<a name="l00475"></a>00475 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size(); i++ )
<a name="l00476"></a>00476 cov[i].resize( <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size(), 0.0 );
<a name="l00477"></a>00477
<a name="l00478"></a>00478 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size(); i++ )
<a name="l00479"></a>00479 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size(); j++ )
<a name="l00480"></a>00480 cov[i][j] = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a68f9abfe3d088353c41bf490c1ef5339" title="The inverse of the quasi-Hessian.">m_M</a>[i][j];
<a name="l00481"></a>00481
<a name="l00482"></a>00482 <span class="comment">// set return flag as EXIT_SUCCESS if cov is Positive Definite,</span>
<a name="l00483"></a>00483 <span class="comment">// and to EXIT_FAILURE otherwise.</span>
<a name="l00484"></a>00484 <span class="keywordtype">int</span> flag = <a class="code" href="namespacehippodraw_1_1Numeric.html#aecbdf19ba7109e37b7e89e490aaa453a" title="The subroutine which does cholesky factorization of a given Symmetric positive definite...">cholFactor</a>( cov );
<a name="l00485"></a>00485
<a name="l00486"></a>00486 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size(); i++ )
<a name="l00487"></a>00487 <span class="keywordflow">for</span>( <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> j = 0; j < <a class="code" href="classhippodraw_1_1BFGSFitter.html#a1b38df992a01ae77f2f06a11c4296feb" title="The initial value to start the iteration from.">m_xinit</a>.size(); j++ )
<a name="l00488"></a>00488 cov[i][j] = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a68f9abfe3d088353c41bf490c1ef5339" title="The inverse of the quasi-Hessian.">m_M</a>[i][j];
<a name="l00489"></a>00489
<a name="l00490"></a>00490 <span class="keywordflow">return</span> flag;
<a name="l00491"></a>00491 }
<a name="l00492"></a>00492
<a name="l00493"></a>00493
<a name="l00494"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a6f92ab21786d5a4f298a2f1dd4782986">00494</a> <span class="keywordtype">double</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a6f92ab21786d5a4f298a2f1dd4782986" title="Given a string, this function returns the value of the associated iteration parameter...">BFGSFitter::iterParam</a> ( std::string name )
<a name="l00495"></a>00495 {
<a name="l00496"></a>00496 <span class="comment">// First check if user is attempting to access the max_iterations</span>
<a name="l00497"></a>00497 <span class="comment">// This is a hack, but it was necessary to ensure a uniform interface</span>
<a name="l00498"></a>00498 <span class="keywordflow">if</span>( name == <span class="stringliteral">"max_iterations"</span> )
<a name="l00499"></a>00499 <span class="keywordflow">return</span> <a class="code" href="classhippodraw_1_1Fitter.html#a35b1e829c1c13a7aef0c1b4cc1895c34" title="The maximum number of iterations allowed in attempting the fit.">m_max_iterations</a>;
<a name="l00500"></a>00500
<a name="l00501"></a>00501 <span class="comment">// Don't use map::operator[]() to find the name and its</span>
<a name="l00502"></a>00502 <span class="comment">// associated value, as it will create one if it doesn't exist.</span>
<a name="l00503"></a>00503 map< string, double * >::const_iterator it
<a name="l00504"></a>00504 = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>.find ( name );
<a name="l00505"></a>00505
<a name="l00506"></a>00506 <span class="keywordflow">if</span> ( it == <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>.end () )
<a name="l00507"></a>00507 cout << name << <span class="stringliteral">" is not a valid iteration parameter name"</span> << endl;
<a name="l00508"></a>00508 <span class="keywordflow">else</span>
<a name="l00509"></a>00509 <span class="keywordflow">return</span> *<a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[name];
<a name="l00510"></a>00510
<a name="l00511"></a>00511 <span class="keywordflow">return</span> 0.0;
<a name="l00512"></a>00512 }
<a name="l00513"></a>00513
<a name="l00514"></a><a class="code" href="classhippodraw_1_1BFGSFitter.html#a09c36569584b6cbfe98aef484581e742">00514</a> <span class="keywordtype">int</span> <a class="code" href="classhippodraw_1_1BFGSFitter.html#a09c36569584b6cbfe98aef484581e742" title="Given a string and a double, this function sets the value of the associated iteration...">BFGSFitter::setIterParam</a> ( std::string name, <span class="keywordtype">double</span> value )
<a name="l00515"></a>00515
<a name="l00516"></a>00516 {
<a name="l00517"></a>00517
<a name="l00518"></a>00518 <span class="comment">// First check if user is attempting to modify the max_iterations</span>
<a name="l00519"></a>00519 <span class="comment">// This is a hack, but it was necessary to ensure a uniform interface</span>
<a name="l00520"></a>00520 <span class="keywordflow">if</span>( name == <span class="stringliteral">"max_iterations"</span> )
<a name="l00521"></a>00521 {
<a name="l00522"></a>00522 <a class="code" href="classhippodraw_1_1Fitter.html#a35b1e829c1c13a7aef0c1b4cc1895c34" title="The maximum number of iterations allowed in attempting the fit.">m_max_iterations</a> = ( int ) value;
<a name="l00523"></a>00523 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00524"></a>00524 }
<a name="l00525"></a>00525
<a name="l00526"></a>00526 <span class="comment">// Now start worrying about the other parameters. </span>
<a name="l00527"></a>00527 <span class="comment">// Don't diretly use map::operator[]() to find the name and its</span>
<a name="l00528"></a>00528 <span class="comment">// associated value, as it will create one if it doesn't exist.</span>
<a name="l00529"></a>00529 map< string, double * >::const_iterator it
<a name="l00530"></a>00530 = <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>.find ( name );
<a name="l00531"></a>00531
<a name="l00532"></a>00532 <span class="keywordflow">if</span> ( it == <a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>.end () )
<a name="l00533"></a>00533 {
<a name="l00534"></a>00534 cout << name << <span class="stringliteral">" is not a valid iteration parameter name"</span> << endl;
<a name="l00535"></a>00535 <span class="keywordflow">return</span> EXIT_FAILURE;
<a name="l00536"></a>00536 }
<a name="l00537"></a>00537 <span class="keywordflow">else</span>
<a name="l00538"></a>00538 {
<a name="l00539"></a>00539 *<a class="code" href="classhippodraw_1_1BFGSFitter.html#a684aa15e32344228d4d248a447a7758a" title="Map of the various iteration parameters to their name.">m_iter_params</a>[name] = value;
<a name="l00540"></a>00540 <span class="keywordflow">return</span> EXIT_SUCCESS;
<a name="l00541"></a>00541 }
<a name="l00542"></a>00542
<a name="l00543"></a>00543 <span class="keywordflow">return</span> EXIT_FAILURE;
<a name="l00544"></a>00544 }
</pre></div></div>
<hr size="1"><address style="align: right;"><small>
Generated for HippoDraw Class Library by <a href="http://www.doxygen.org/index.html">
<img src="doxygen.png" alt="doxygen" align="middle" border=0>
</a> </small></address>
</body>
</html>
