nojhan/eodev
Archived
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

code formating

Browse source
This commit is contained in:
nojhan 2011-09-15 17:27:46 +02:00
commit 37a2c68b69
4 changed files with 142 additions and 192 deletions
Download patch file Download diff file Expand all files Collapse all files

160
edo/src/edoSamplerNormalMulti.h
View file

@ -40,136 +40,138 @@ class edoSamplerNormalMulti : public edoSampler< edoNormalMulti< EOT > >
public:
typedef typename EOT::AtomType AtomType;
edoSamplerNormalMulti( edoRepairer<EOT> & repairer ) : edoSampler( repairer) {}
class Cholesky
{
public:
Cholesky( const ublas::symmetric_matrix< AtomType, ublas::lower >& V)
{
unsigned int Vl = V.size1();
Cholesky( const ublas::symmetric_matrix< AtomType, ublas::lower >& V)
{
unsigned int Vl = V.size1();
assert(Vl > 0);
assert(Vl > 0);
unsigned int Vc = V.size2();
unsigned int Vc = V.size2();
assert(Vc > 0);
assert(Vc > 0);
assert( Vl == Vc );
assert( Vl == Vc );
_L.resize(Vl);
_L.resize(Vl);
unsigned int i,j,k;
unsigned int i,j,k;
// first column
i=0;
// first column
i=0;
// diagonal
j=0;
_L(0, 0) = sqrt( V(0, 0) );
// diagonal
j=0;
_L(0, 0) = sqrt( V(0, 0) );
// end of the column
for ( j = 1; j < Vc; ++j )
{
_L(j, 0) = V(0, j) / _L(0, 0);
}
// end of the column
for ( j = 1; j < Vc; ++j )
{
_L(j, 0) = V(0, j) / _L(0, 0);
}
// end of the matrix
for ( i = 1; i < Vl; ++i ) // each column
{
// end of the matrix
for ( i = 1; i < Vl; ++i ) // each column
{
// diagonal
double sum = 0.0;
// diagonal
double sum = 0.0;
for ( k = 0; k < i; ++k)
{
sum += _L(i, k) * _L(i, k);
}
for ( k = 0; k < i; ++k)
{
sum += _L(i, k) * _L(i, k);
}
_L(i,i) = sqrt( fabs( V(i,i) - sum) );
_L(i,i) = sqrt( fabs( V(i,i) - sum) );
for ( j = i + 1; j < Vl; ++j ) // rows
{
// one element
sum = 0.0;
for ( j = i + 1; j < Vl; ++j ) // rows
{
// one element
sum = 0.0;
for ( k = 0; k < i; ++k )
{
sum += _L(j, k) * _L(i, k);
}
for ( k = 0; k < i; ++k )
{
sum += _L(j, k) * _L(i, k);
}
_L(j, i) = (V(j, i) - sum) / _L(i, i);
}
}
}
_L(j, i) = (V(j, i) - sum) / _L(i, i);
}
}
}
const ublas::symmetric_matrix< AtomType, ublas::lower >& get_L() const {return _L;}
const ublas::symmetric_matrix< AtomType, ublas::lower >& get_L() const {return _L;}
private:
ublas::symmetric_matrix< AtomType, ublas::lower > _L;
ublas::symmetric_matrix< AtomType, ublas::lower > _L;
};
edoSamplerNormalMulti( edoBounder< EOT > & bounder )
: edoSampler< edoNormalMulti< EOT > >( bounder )
: edoSampler< edoNormalMulti< EOT > >( bounder )
{}
EOT sample( edoNormalMulti< EOT >& distrib )
{
unsigned int size = distrib.size();
unsigned int size = distrib.size();
assert(size > 0);
assert(size > 0);
//-------------------------------------------------------------
// Cholesky factorisation gererating matrix L from covariance
// matrix V.
// We must use cholesky.get_L() to get the resulting matrix.
//
// L = cholesky decomposition of varcovar
//-------------------------------------------------------------
//-------------------------------------------------------------
// Cholesky factorisation gererating matrix L from covariance
// matrix V.
// We must use cholesky.get_L() to get the resulting matrix.
//
// L = cholesky decomposition of varcovar
//-------------------------------------------------------------
Cholesky cholesky( distrib.varcovar() );
ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L();
Cholesky cholesky( distrib.varcovar() );
ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L();
//-------------------------------------------------------------
//-------------------------------------------------------------
//-------------------------------------------------------------
// T = vector of size elements drawn in N(0,1) rng.normal(1.0)
//-------------------------------------------------------------
//-------------------------------------------------------------
// T = vector of size elements drawn in N(0,1) rng.normal(1.0)
//-------------------------------------------------------------
ublas::vector< AtomType > T( size );
ublas::vector< AtomType > T( size );
for ( unsigned int i = 0; i < size; ++i )
{
T( i ) = rng.normal( 1.0 );
}
for ( unsigned int i = 0; i < size; ++i )
{
T( i ) = rng.normal( 1.0 );
}
//-------------------------------------------------------------
//-------------------------------------------------------------
//-------------------------------------------------------------
// LT = prod( L, T )
//-------------------------------------------------------------
//-------------------------------------------------------------
// LT = prod( L, T )
//-------------------------------------------------------------
ublas::vector< AtomType > LT = ublas::prod( L, T );
ublas::vector< AtomType > LT = ublas::prod( L, T );
//-------------------------------------------------------------
//-------------------------------------------------------------
//-------------------------------------------------------------
// solution = means + LT
//-------------------------------------------------------------
//-------------------------------------------------------------
// solution = means + LT
//-------------------------------------------------------------
ublas::vector< AtomType > mean = distrib.mean();
ublas::vector< AtomType > mean = distrib.mean();
ublas::vector< AtomType > ublas_solution = mean + LT;
ublas::vector< AtomType > ublas_solution = mean + LT;
EOT solution( size );
EOT solution( size );
std::copy( ublas_solution.begin(), ublas_solution.end(), solution.begin() );
std::copy( ublas_solution.begin(), ublas_solution.end(), solution.begin() );
//-------------------------------------------------------------
//-------------------------------------------------------------
return solution;
return solution;
}
};