|
Sleipnir
|
00001 /***************************************************************************** 00002 * This file is provided under the Creative Commons Attribution 3.0 license. 00003 * 00004 * You are free to share, copy, distribute, transmit, or adapt this work 00005 * PROVIDED THAT you attribute the work to the authors listed below. 00006 * For more information, please see the following web page: 00007 * http://creativecommons.org/licenses/by/3.0/ 00008 * 00009 * This file is a component of the Sleipnir library for functional genomics, 00010 * authored by: 00011 * Curtis Huttenhower (chuttenh@princeton.edu) 00012 * Mark Schroeder 00013 * Maria D. Chikina 00014 * Olga G. Troyanskaya (ogt@princeton.edu, primary contact) 00015 * 00016 * If you use this library, the included executable tools, or any related 00017 * code in your work, please cite the following publication: 00018 * Curtis Huttenhower, Mark Schroeder, Maria D. Chikina, and 00019 * Olga G. Troyanskaya. 00020 * "The Sleipnir library for computational functional genomics" 00021 *****************************************************************************/ 00022 #include "stdafx.h" 00023 #include "bayesnet.h" 00024 #include "dataset.h" 00025 00026 #ifndef NO_SMILE 00027 00028 namespace Sleipnir { 00029 00030 bool CBayesNetSmileImpl::LearnELR( const IDataset* pData, size_t iIterations, bool fZero ) { 00031 size_t i, j, iParameters; 00032 TVecVecD vecvecfBeta, vecvecfGradient, vecvecfPrev, vecvecfDirection, vecvecfOriginal; 00033 DSL_Dmatrix* pMatrix; 00034 DSL_intArray veciCoords; 00035 float dAX, dBX, dDotNew, dDotOld, dDotON, dB; 00036 // TTrieData TrieData( pData, 0 ); 00037 TMapData mapData; 00038 vector<bool> vecfHidden; 00039 00040 if( !m_fSmileNet ) 00041 return false; 00042 00043 vecfHidden.resize( pData->GetExperiments( ) ); 00044 for( i = 0; i < pData->GetExperiments( ); ++i ) 00045 vecfHidden[ i ] = pData->IsHidden( i ); 00046 EncodeData( pData, mapData ); 00047 iParameters = ELRCountParameters( ); 00048 vecvecfBeta.resize( m_SmileNet.GetNumberOfNodes( ) ); 00049 vecvecfGradient.resize( m_SmileNet.GetNumberOfNodes( ) ); 00050 vecvecfPrev.resize( m_SmileNet.GetNumberOfNodes( ) ); 00051 vecvecfDirection.resize( m_SmileNet.GetNumberOfNodes( ) ); 00052 vecvecfOriginal.resize( m_SmileNet.GetNumberOfNodes( ) ); 00053 for( i = 0; i < vecvecfBeta.size( ); ++i ) { 00054 pMatrix = m_SmileNet.GetNode( (int)i )->Definition( )->GetMatrix( ); 00055 vecvecfOriginal[ i ].resize( pMatrix->GetSize( ) ); 00056 vecvecfDirection[ i ].resize( pMatrix->GetSize( ) ); 00057 vecvecfPrev[ i ].resize( pMatrix->GetSize( ) ); 00058 vecvecfGradient[ i ].resize( pMatrix->GetSize( ) ); 00059 vecvecfBeta[ i ].resize( pMatrix->GetSize( ) ); 00060 for( j = 0; j < vecvecfBeta[ i ].size( ); ++j ) 00061 vecvecfBeta[ i ][ j ] = (float)log( (*pMatrix)[ (int)j ] ); } 00062 dAX = 0; 00063 dBX = 0.1f; 00064 00065 copy( vecvecfBeta.begin( ), vecvecfBeta.end( ), vecvecfOriginal.begin( ) ); 00066 ELRComputeGradient( vecfHidden, mapData, fZero, vecvecfGradient ); 00067 copy( vecvecfGradient.begin( ), vecvecfGradient.end( ), vecvecfPrev.begin( ) ); 00068 copy( vecvecfGradient.begin( ), vecvecfGradient.end( ), vecvecfDirection.begin( ) ); 00069 dDotNew = ELRDot( vecvecfGradient, vecvecfPrev ); 00070 00071 for( i = 0; i < iIterations; ++i ) { 00072 g_CatSleipnir( ).notice( "CBayesNetSmile::LearnELR( %d, %d ) iteration %d/%d", 00073 iIterations, fZero, i, iIterations ); 00074 if( !dDotNew ) 00075 continue; 00076 ELRNormalizeDirection( vecvecfDirection ); 00077 ELRLineSearch( vecfHidden, mapData, vecvecfDirection, vecvecfOriginal, vecvecfBeta, dAX, dBX, fZero ); 00078 copy( vecvecfBeta.begin( ), vecvecfBeta.end( ), vecvecfOriginal.begin( ) ); 00079 ELRCopyParameters( vecvecfBeta ); 00080 ELRComputeGradient( vecfHidden, mapData, fZero, vecvecfGradient ); 00081 00082 dDotOld = dDotNew; 00083 dDotON = ELRDot( vecvecfGradient, vecvecfPrev ); 00084 copy( vecvecfGradient.begin( ), vecvecfGradient.end( ), vecvecfPrev.begin( ) ); 00085 dDotNew = ELRDot( vecvecfGradient, vecvecfPrev ); 00086 dB = ( dDotNew - dDotON ) / dDotOld; 00087 if( !( ( i + 1 ) % iParameters ) || ( dB <= 0 ) ) 00088 copy( vecvecfGradient.begin( ), vecvecfGradient.end( ), vecvecfDirection.begin( ) ); 00089 else 00090 ELRComputeDirection( dB, vecvecfGradient, vecvecfDirection ); } 00091 00092 return true; } 00093 00094 size_t CBayesNetSmileImpl::ELRCountParameters( ) const { 00095 int i; 00096 size_t iRet; 00097 00098 for( iRet = i = 0; i < m_SmileNet.GetNumberOfNodes( ); ++i ) 00099 iRet += m_SmileNet.GetNode( i )->Definition( )->GetSize( ); 00100 00101 return iRet; } 00102 00103 void CBayesNetSmileImpl::ELRCopyParameters( TVecVecD& vecvecfBeta ) { 00104 static const float c_dMinimum = log( FLT_MIN ); 00105 size_t i, j, k, iIndex, iDomains; 00106 float dSum, dCur, dMax; 00107 DSL_nodeDefinition* pDef; 00108 DSL_Dmatrix* pDefs; 00109 00110 for( i = 0; i < vecvecfBeta.size( ); ++i ) { 00111 pDef = m_SmileNet.GetNode( (int)i )->Definition( ); 00112 pDefs = pDef->GetMatrix( ); 00113 iDomains = pDef->GetNumberOfOutcomes( ); 00114 for( j = 0; j < vecvecfBeta[ i ].size( ); j += iDomains ) { 00115 dSum = 0; 00116 dMax = vecvecfBeta[ i ][ j ]; 00117 for( k = 1; k < iDomains; ++k ) 00118 if( ( dCur = vecvecfBeta[ i ][ j + k ] ) > dMax ) 00119 dMax = dCur; 00120 for( k = 0; k < iDomains; ++k ) { 00121 iIndex = j + k; 00122 if( ( vecvecfBeta[ i ][ iIndex ] -= dMax ) < c_dMinimum ) 00123 vecvecfBeta[ i ][ iIndex ] = c_dMinimum; 00124 dSum += exp( vecvecfBeta[ i ][ iIndex ] ); } 00125 if( dSum ) 00126 for( k = 0; k < iDomains; ++k ) { 00127 iIndex = j + k; 00128 (*pDefs)[ (int)iIndex ] = exp( vecvecfBeta[ i ][ iIndex ] ) / dSum; } } } } 00129 00130 // EXPENSIVE 00131 void CBayesNetSmileImpl::ELRComputeGradient( const vector<bool>& vecfHidden, const TMapData& mapData, bool fZero, 00132 TVecVecD& vecvecfGradient ) { 00133 size_t i, j; 00134 00135 for( i = 0; i < vecvecfGradient.size( ); ++i ) 00136 for( j = 0; j < vecvecfGradient[ i ].size( ); ++j ) 00137 vecvecfGradient[ i ][ j ] = 0; 00138 00139 for( TMapData::const_iterator iterData = mapData.begin( ); iterData != mapData.end( ); ++iterData ) 00140 if( IsAnswer( iterData->first ) && FillCPTs( vecfHidden, iterData->first, fZero, false ) ) { 00141 ELRUpdateGradient( -(float)iterData->second, vecvecfGradient ); 00142 m_SmileNet.GetNode( 0 )->Value( )->SetEvidence( iterData->first[ 0 ] - c_cBase ); 00143 ELRUpdateGradient( (float)iterData->second, vecvecfGradient ); } } 00144 /* 00145 for( TTrieData::iterator IterData( TrieData ); !IterData.IsDone( ); IterData.Next( ) ) 00146 if( IterData.GetPosition( )[ 0 ] && FillCPTs( vecfHidden, IterData.GetPosition( ), fZero, false ) ) { 00147 ELRUpdateGradient( -(float)IterData.Get( ), vecvecfGradient ); 00148 m_SmileNet.GetNode( 0 )->Value( )->SetEvidence( IterData.GetPosition( )[ 0 ] - 1 ); 00149 ELRUpdateGradient( (float)IterData.Get( ), vecvecfGradient ); } } 00150 */ 00151 00152 void CBayesNetSmileImpl::ELRUpdateGradient( float dRate, TVecVecD& vecvecfGradient ) { 00153 size_t i, j; 00154 double dSum, dPosterior; 00155 int k, iTmp, iDomain, iDomains, iEvidence; 00156 DSL_node* pNode; 00157 DSL_nodeValue* pValue; 00158 DSL_nodeDefinition* pDef; 00159 DSL_Dmatrix* pValues; 00160 DSL_Dmatrix* pDefs; 00161 DSL_intArray veciParents, veciCoords; 00162 00163 m_SmileNet.UpdateBeliefs( ); 00164 for( i = 0; i < vecvecfGradient.size( ); ++i ) { 00165 pNode = m_SmileNet.GetNode( (int)i ); 00166 veciParents = pNode->Parents( ); 00167 pValue = pNode->Value( ); 00168 iEvidence = pValue->GetEvidence( ); 00169 pValues = pValue->GetMatrix( ); 00170 pDef = pNode->Definition( ); 00171 iDomains = pDef->GetNumberOfOutcomes( ); 00172 pDefs = pDef->GetMatrix( ); 00173 00174 dSum = 0; 00175 for( j = 0,pDefs->IndexToCoordinates( iTmp = 0, veciCoords ); iTmp != DSL_OUT_OF_RANGE; 00176 ++j,iTmp = pDefs->NextCoordinates( veciCoords ) ) { 00177 iDomain = veciCoords[ veciCoords.GetSize( ) - 1 ]; 00178 if( veciParents.NumItems( ) ) { 00179 if( iEvidence == DSL_OUT_OF_RANGE ) 00180 dPosterior = (*pValues)[ iDomain ]; 00181 else if( iDomain == iEvidence ) { 00182 dPosterior = 1; 00183 for( k = 0; k < veciParents.NumItems( ); ++k ) 00184 dPosterior *= m_SmileNet.GetNode( veciParents[ k ] )->Value( 00185 )->GetMatrix( )->Subscript( veciCoords[ k ] ); } 00186 else 00187 dPosterior = 0; } 00188 else 00189 dPosterior = (*pValues)[ veciCoords[ 0 ] ]; 00190 dSum += dPosterior; 00191 vecvecfGradient[ i ][ j ] += dRate * (float)dPosterior; 00192 00193 if( ( iDomain + 1 ) == iDomains ) { 00194 for( k = 0; k < iDomains; ++k ) { 00195 veciCoords[ veciCoords.GetSize( ) - 1 ] = k; 00196 vecvecfGradient[ i ][ k + j - iDomains + 1 ] -= 00197 (float)( dRate * dSum * (*pDefs)[ veciCoords ] ); } 00198 dSum = 0; } } } } 00199 00200 float CBayesNetSmileImpl::ELRDot( const TVecVecD& vecvecfOne, const TVecVecD& vecvecfTwo ) { 00201 size_t i, j; 00202 float dRet; 00203 00204 dRet = 0; 00205 for( i = 0; i < vecvecfOne.size( ); ++i ) 00206 for( j = 0; j < vecvecfOne[ i ].size( ); ++j ) 00207 dRet += vecvecfOne[ i ][ j ] * vecvecfTwo[ i ][ j ]; 00208 00209 return dRet; } 00210 00211 void CBayesNetSmileImpl::ELRNormalizeDirection( TVecVecD& vecvecfDirection ) const { 00212 size_t i, j; 00213 float d, dSum; 00214 int k, iDomains; 00215 00216 for( i = 0; i < vecvecfDirection.size( ); ++i ) { 00217 iDomains = m_SmileNet.GetNode( (int)i )->Definition( )->GetNumberOfOutcomes( ); 00218 for( j = 0; j < vecvecfDirection[ i ].size( ); j += iDomains ) { 00219 dSum = 0; 00220 for( k = 0; k < iDomains; ++k ) { 00221 d = vecvecfDirection[ i ][ j + k ]; 00222 dSum += d * d; } 00223 00224 if( dSum = sqrt( dSum ) ) 00225 for( k = 0; k < iDomains; ++k ) 00226 vecvecfDirection[ i ][ j + k ] /= dSum; } } } 00227 00228 // EXPENSIVE 00229 float CBayesNetSmileImpl::ELRLineSearch( const vector<bool>& vecfHidden, const TMapData& mapData, 00230 const TVecVecD& vecvecfDirection, const TVecVecD& vecvecfOriginal, TVecVecD& vecvecfBeta, 00231 float& dAX, float& dBX, bool fZero ) { 00232 float dFA, dFB, dFC, dCX; 00233 00234 ELRBracket( vecfHidden, mapData, vecvecfDirection, vecvecfOriginal, vecvecfBeta, 00235 dAX, dBX, dCX, dFA, dFB, dFC, fZero ); 00236 return ELRBrent( vecfHidden, mapData, vecvecfDirection, vecvecfOriginal, vecvecfBeta, 00237 dAX, dBX, dCX, dFA, dFB, dFC, fZero ); } 00238 00239 float CBayesNetSmileImpl::ELREvalFunction( const vector<bool>& vecfHidden, const TMapData& mapData, float dX, 00240 const TVecVecD& vecvecfDirection, const TVecVecD& vecvecfOriginal, TVecVecD& vecvecfBeta, bool fZero ) { 00241 size_t i, j; 00242 00243 for( i = 0; i < vecvecfBeta.size( ); ++i ) 00244 for( j = 0; j < vecvecfBeta[ i ].size( ); ++j ) 00245 vecvecfBeta[ i ][ j ] = ( dX * vecvecfDirection[ i ][ j ] ) + vecvecfOriginal[ i ][ j ]; 00246 ELRCopyParameters( vecvecfBeta ); 00247 00248 return -ELRConditionalLikelihood( vecfHidden, mapData, fZero ); } 00249 00250 float CBayesNetSmileImpl::ELRAvoidZero( float d ) { 00251 static const float c_dTiny = 1e-10f; 00252 00253 if( d >= 0 ) { 00254 if( d < c_dTiny ) 00255 return c_dTiny; } 00256 else if( d > -c_dTiny ) 00257 return -c_dTiny; 00258 00259 return d; } 00260 00261 void CBayesNetSmileImpl::ELRBracket( const vector<bool>& vecfHidden, const TMapData& mapData, 00262 const TVecVecD& vecvecfDirection, const TVecVecD& vecvecfOriginal, TVecVecD& vecvecfBeta, float& dAX, 00263 float& dBX, float& dCX, float& dFA, float& dFB, float& dFC, bool fZero ) { 00264 static const float c_dGolden = 1.618034f; 00265 static const float c_dLimit = 100; 00266 static const size_t c_iIterations = 100; 00267 size_t i; 00268 float dR, dQ, dU, dFU, dULim; 00269 00270 dFA = ELREvalFunction( vecfHidden, mapData, dAX, vecvecfDirection, vecvecfOriginal, vecvecfBeta, fZero ); 00271 dFB = ELREvalFunction( vecfHidden, mapData, dBX, vecvecfDirection, vecvecfOriginal, vecvecfBeta, fZero ); 00272 if( dFB > dFA ) { 00273 swap( dAX, dBX ); 00274 swap( dFA, dFB ); } 00275 00276 dCX = dBX + ( c_dGolden * ( dBX - dAX ) ); 00277 dFC = ELREvalFunction( vecfHidden, mapData, dCX, vecvecfDirection, vecvecfOriginal, vecvecfBeta, fZero ); 00278 for( i = 0; ( i < c_iIterations ) && ( dFB > dFC ); ++i ) { 00279 /* finding the minimum point dU of the extrapolated parabola*/ 00280 dR = ( dBX - dAX ) * ( dFB - dFC ); 00281 dQ = ( dBX - dCX ) * ( dFB - dFA ); 00282 dU = dBX - ( ( ( ( dBX - dCX ) * dQ ) - ( ( dBX - dAX ) * dR ) ) / ( 2 * ELRAvoidZero( dQ - dR ) ) ); 00283 dULim = dBX + ( c_dLimit * ( dCX - dBX ) ); 00284 00285 /* bx>u and u>cx, or u is between b and c */ 00286 if( ( ( dBX - dU ) * ( dU - dCX ) ) > 0 ) { 00287 dFU = ELREvalFunction( vecfHidden, mapData, dU, vecvecfDirection, vecvecfOriginal, vecvecfBeta, fZero ); 00288 00289 /* the minimum is between b and c, because c is lower than b, 00290 * otherwise we can not be in the loop. So u is lower than both b 00291 * and c. We are done. */ 00292 if( dFU < dFC ) { 00293 dAX = dBX; 00294 dFA = dFB; 00295 dBX = dU; 00296 dFB = dFU; 00297 break; } 00298 /* b is lower than both a and u. we are done*/ 00299 else if( dFU > dFB ) { 00300 dCX = dU; 00301 dFC = dFU; 00302 break; } 00303 /* u is lower than b, but higher than c. useless. 00304 * extend it beyound c to prepare for the next iteration 00305 * u will be copied to c. */ 00306 00307 dU = dCX + ( c_dGolden * ( dCX - dBX ) ); 00308 dFU = ELREvalFunction( vecfHidden, mapData, dU, vecvecfDirection, vecvecfOriginal, vecvecfBeta, fZero ); } 00309 /* u is between c and the far limit */ 00310 else if( ( ( dCX - dU ) * ( dU - dULim ) ) > 0 ) { 00311 dFU = ELREvalFunction( vecfHidden, mapData, dU, vecvecfDirection, vecvecfOriginal, vecvecfBeta, fZero ); 00312 /* if u is still lower than c, extend the points beyond u and shift the values */ 00313 if( dFU < dFC ) { 00314 dBX = dCX; 00315 dCX = dU; 00316 dU = dCX + ( c_dGolden * ( dCX - dBX ) ); 00317 dFB = dFC; 00318 dFC = dFU; 00319 dFU = ELREvalFunction( vecfHidden, mapData, dU, vecvecfDirection, vecvecfOriginal, vecvecfBeta, 00320 fZero ); } } 00321 /* u is beyond the far limit. ulim is between u and c. c and ulim defines the direction. */ 00322 else if( ( ( dU - dULim ) * ( dULim - dCX ) ) > 0 ) 00323 dFU = ELREvalFunction( vecfHidden, mapData, dU = dULim, vecvecfDirection, vecvecfOriginal, vecvecfBeta, 00324 fZero ); 00325 /* other cases... what could be left? */ 00326 else 00327 dFU = ELREvalFunction( vecfHidden, mapData, dU = dCX + ( c_dGolden * ( dCX - dBX ) ), vecvecfDirection, 00328 vecvecfOriginal, vecvecfBeta, fZero ); 00329 00330 /*eliminate oldest point and repeat the process*/ 00331 dAX = dBX; 00332 dBX = dCX; 00333 dCX = dU; 00334 dFA = dFB; 00335 dFB = dFC; 00336 dFC = dFU; } } 00337 00338 float CBayesNetSmileImpl::ELRBrent( const vector<bool>& vecfHidden, const TMapData& mapData, 00339 const TVecVecD& vecvecfDirection, const TVecVecD& vecvecfOriginal, TVecVecD& vecvecfBeta, float& dAX, 00340 float& dBX, float dCX, float dFA, float dFB, float dFC, bool fZero ) { 00341 static const size_t c_iIterations = 100; 00342 static const float c_dZEPS = 1e-8f; 00343 static const float c_dIGolden = 0.3819660f; 00344 static const float c_dTolerance = 1e-4f; 00345 /* d, e, etemp are the step. e is the step before d, etemp is the step 00346 * before e 00347 * tol,tol1,tol2 are the error tolerance. for simple thinking, consider them 0. */ 00348 float a,b,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u=0,v,w,x,xm,e,d; 00349 size_t iter; 00350 00351 /* a, b brackets the minimum. ensure a<=b*/ 00352 if( dAX < dCX ) { 00353 a = dAX; 00354 b = dCX; } 00355 else { 00356 a = dCX; 00357 b = dAX; } 00358 00359 /* initialize all three lowest points to the middle variable of the bracket*/ 00360 e = d = b - a; 00361 x = dBX; 00362 fx = dFB; 00363 if( dFA < dFC ) { 00364 w = dAX; 00365 fw = dFA; 00366 v = dCX; 00367 fv = dFC; } 00368 else { 00369 w = dCX; 00370 fw = dFC; 00371 v = dAX; 00372 fv = dFA; } 00373 00374 for( iter = 0; iter < c_iIterations; ++iter ) { 00375 xm = 0.5f * ( a + b ); 00376 /* error tolerance:2*tol*x. note that tol, tol1, tol2 are all positive*/ 00377 tol2 = 2 * ( tol1 = ( ( c_dTolerance * fabs( x ) ) + c_dZEPS ) ); 00378 /* the sum of x - xm and 0.5*(b-a) smaller than tolerance*/ 00379 if( ( fw - fx ) <= ( c_dTolerance * fx ) ) { 00380 dAX = a; 00381 dBX = b; 00382 return fx; } 00383 00384 /* e is used to check to see whether it is ok to take a parabolic step 00385 * parabolic step will not be considered if the step before previous 00386 * step was too small , made no progress. 00387 * I'm not sure when could this happen... if progress is too small, 00388 * it should be done.... unless it's the first time */ 00389 if( fabs( e ) > tol1 ) { 00390 /* fit a parabola*/ 00391 r = ( x - w ) * ( fx - fv ); 00392 q = ( x - v ) * ( fx - fw ); 00393 p = ( ( x - v ) * q ) - ( ( x - w ) * r ); 00394 q = 2 * ( q - r ); 00395 00396 /*ensure q is positive, and reverse the ratio p/q 00397 * these extra steps are used to avoid dividing by zero, making the 00398 * denominator positive ensure the later inequality holds because it is 00399 * multiplied by a positive number */ 00400 if( q > 0 ) 00401 p *= -1; 00402 else 00403 q = -q; 00404 00405 etemp = e; 00406 /*the previouse step size*/ 00407 e = d; 00408 00409 /* after this iteration, e will be either the previous step size 00410 * (if the current step is a parabolic step), or the size of the 00411 * current larger segment, (if the current step is a GOLD step). 00412 * therefore, on the next iteration, a parabolic step will be 00413 * taken only if the parabolic step is smaller than the 00414 * previous step, or smaller than the current larger segment */ 00415 00416 /* the first part test if the current step is greater than the step 00417 * before the previous step. If yes, it may not be making progress 00418 * and may not converge using parabola, so use GOLD step instead. 00419 * 00420 * the 2nd part test the new minimum of the parabola x+p/q (which 00421 * is the old x-p/q) is outside of (a,b) 00422 * e is set to be the larger of the 2 segments */ 00423 if( ( fabs( p ) >= fabs( 0.5 * q * etemp ) ) || 00424 ( ( p <= ( q * ( a - x ) ) ) || ( p >= ( q * ( b - x ) ) ) ) ) 00425 d = c_dIGolden * ( e = ( x >= xm ? ( a - x ) : ( b - x ) ) ); 00426 /* parabola is good. update the minimum u */ 00427 else { 00428 /* parabolic step*/ 00429 d = p / q; 00430 /* notice that u's value is not fixed yet, it's only set for 00431 * testing purpose now. it will be fixed later 00432 * Bad coding. Here should really make the code cleaner by using another 00433 * variable, say uTest, or boundaryTest. */ 00434 u = x + d; 00435 00436 /* if the minimum is too close to a or b, step away by tol1 to 00437 * ensure there is progress made.*/ 00438 if( ( u - a ) < tol2 ) 00439 d += tol1; 00440 else if( ( b - u ) < tol2 ) 00441 d -= tol1; } } 00442 /* take a GOLD step in case of problem*/ 00443 else 00444 d = c_dIGolden * ( e = ( ( x > xm ) ? ( a - x ) : ( b - x ) ) ); 00445 00446 /* update the new minimum 00447 * if d<tolerance, step by tolerance, because otherwise it would give 00448 * the "same" result and make no progress */ 00449 if( fabs( d ) >=tol1 ) 00450 u = x + d; 00451 else { 00452 if( d > 0 ) 00453 u = x + tol1; 00454 else 00455 u = x - tol1; } 00456 fu = ELREvalFunction( vecfHidden, mapData, u, vecvecfDirection, vecvecfOriginal, vecvecfBeta, fZero ); 00457 00458 /* u is the lowest point of the function, 00459 * lower than the previous lowest point x. 00460 * so make x one of the boundary points. 00461 * then old point v will be discarded. 00462 * and one of old a,b will be discarded, and will be the 2nd lowest 00463 * point w, which was the old lowest point x. 00464 * so w will be one of the boundary points a,b. */ 00465 if( fu < fx ) { 00466 if( u >= x ) 00467 a = x; 00468 else 00469 b = x; 00470 /*update three lowest points x,v,w*/ 00471 v = w; 00472 w = x; 00473 x = u; 00474 fv = fw; 00475 fw = fx; 00476 fx = fu; } 00477 /* x is still the lowest point*/ 00478 else { 00479 /* we know u is between a and b from before, 00480 * so we can reduce the bracket */ 00481 if( u < x ) 00482 a = u; 00483 else 00484 b = u; 00485 00486 /* if u is lower than 2nd lowest point w, make u the 2nd lowest 00487 * point, and w the 3rd 00488 * if the lowest point is same as the 2nd lowest point */ 00489 if( ( fu <= fw ) || ( w == x ) ) { 00490 v = w; 00491 fv = fw; 00492 w = u; 00493 fw = fu; } 00494 /* if u is lower than the 3rd lowest point v, 00495 * or if the 3rd lowest point is same as the lowest or 2nd lowest point */ 00496 else if( ( fu <= fv ) || ( v == x ) || ( v == w ) ) { 00497 v = u; 00498 fv = fu; } } } 00499 00500 g_CatSleipnir( ).warn( "CBayesNetSmileImpl::ELRBrent( ) too many iterations" ); 00501 dAX = a; 00502 dBX = b; 00503 return fx; } 00504 00505 float CBayesNetSmileImpl::ELRConditionalLikelihood( const vector<bool>& vecfHidden, const TMapData& mapData, 00506 bool fZero ) { 00507 size_t iCount; 00508 float dRet; 00509 DSL_Dmatrix* pMatrix; 00510 00511 iCount = 0; 00512 dRet = 0; 00513 for( TMapData::const_iterator iterData = mapData.begin( ); iterData != mapData.end( ); ++iterData ) 00514 if( IsAnswer( iterData->first ) && FillCPTs( vecfHidden, iterData->first, fZero, false ) ) { 00515 iCount += iterData->second; 00516 m_SmileNet.UpdateBeliefs( ); 00517 pMatrix = m_SmileNet.GetNode( 0 )->Value( )->GetMatrix( ); 00518 dRet += (float)( iterData->second * log( (*pMatrix)[ iterData->first[ 0 ] - c_cBase ] ) ); } 00519 /* 00520 for( TTrieData::iterator IterData( TrieData ); !IterData.IsDone( ); IterData.Next( ) ) 00521 if( IterData.GetPosition( )[ 0 ] && FillCPTs( vecfHidden, IterData.GetPosition( ), fZero, false ) ) { 00522 iCount += IterData.Get( ); 00523 m_SmileNet.UpdateBeliefs( ); 00524 pMatrix = m_SmileNet.GetNode( 0 )->Value( )->GetMatrix( ); 00525 dRet += (float)( IterData.Get( ) * log( (*pMatrix)[ IterData.GetPosition( )[ 0 ] - 1 ] ) ); } 00526 */ 00527 00528 return ( dRet / iCount ); } 00529 00530 void CBayesNetSmileImpl::ELRComputeDirection( float dB, const TVecVecD& vecvecfGradient, 00531 TVecVecD& vecvecfDirection ) { 00532 size_t i, j; 00533 00534 for( i = 0; i < vecvecfDirection.size( ); ++i ) 00535 for( j = 0; j < vecvecfDirection[ i ].size( ); ++j ) 00536 vecvecfDirection[ i ][ j ] = vecvecfGradient[ i ][ j ] + ( dB * vecvecfDirection[ i ][ j ] ); } 00537 00538 } 00539 00540 #endif // NO_SMILE
1.7.6.1