an intermediate commit to keep updates

src/doSamplerNormalMulti.h

@@ -0,0 +1,160 @@
+#ifndef _doSamplerNormalMulti_h
+#define _doSamplerNormalMulti_h
+
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+
+#include <utils/eoRNG.h>
+
+#include "doSampler.h"
+#include "doNormalMulti.h"
+#include "doBounder.h"
+#include "doStats.h"
+
+/**
+ * doSamplerNormalMulti
+ * This class uses the Normal distribution parameters (bounds) to return
+ * a random position used for population sampling.
+ */
+template < typename EOT >
+class doSamplerNormalMulti : public doSampler< doNormalMulti< EOT > >
+{
+public:
+    typedef typename EOT::AtomType AtomType;
+
+    class Cholesky
+    {
+    public:
+	virtual void update( const ublas::symmetric_matrix< AtomType, ublas::lower >& V)
+	{
+	    unsigned int Vl = V.size1();
+
+	    assert(Vl > 0);
+
+	    unsigned int Vc = V.size2();
+
+	    assert(Vc > 0);
+
+	    _L.resize(Vl, Vc);
+
+	    unsigned int i,j,k;
+
+	    // first column
+	    i=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 matrix
+	    for ( i = 1; i < Vl; ++i ) // each column
+		{
+
+		    // diagonal
+		    double sum = 0.0;
+
+		    for ( k = 0; k < i; ++k)
+			{
+			    sum += _L(i, k) * _L(i, k);
+			}
+
+		    assert( ( V(i, i) - sum ) > 0 );
+
+		    _L(i, i) = sqrt( V(i, i) - sum );
+
+		    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);
+				}
+
+			    _L(j, i) = (V(j, i) - sum) / _L(i, i);
+			}
+		}
+	}
+
+	const ublas::symmetric_matrix< AtomType, ublas::lower >& get_L() const {return _L;}
+
+    private:
+	ublas::symmetric_matrix< AtomType, ublas::lower > _L;
+    };
+
+    doSamplerNormalMulti( doBounder< EOT > & bounder )
+	: doSampler< doNormalMulti< EOT > >( bounder )
+    {}
+
+    EOT sample( doNormalMulti< EOT >& distrib )
+    {
+	unsigned int size = distrib.size();
+
+	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 cholesky;
+	cholesky.update( 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)
+	//-------------------------------------------------------------
+
+	ublas::vector< AtomType > T( size );
+
+	for ( unsigned int i = 0; i < size; ++i )
+	    {
+		T( i ) = rng.normal( 1.0 );
+	    }
+
+	//-------------------------------------------------------------
+
+
+	//-------------------------------------------------------------
+	// LT = prod( L, T )
+	//-------------------------------------------------------------
+
+	ublas::vector< AtomType > LT = ublas::prod( L, T );
+
+	//-------------------------------------------------------------
+
+
+	//-------------------------------------------------------------
+	// solution = means + LT
+	//-------------------------------------------------------------
+
+	ublas::vector< AtomType > mean = distrib.mean();
+
+	ublas::vector< AtomType > ublas_solution = mean + LT;
+
+	EOT solution( size );
+
+	std::copy( ublas_solution.begin(), ublas_solution.end(), solution.begin() );
+
+	//-------------------------------------------------------------
+
+	return solution;
+    }
+};
+
+#endif // !_doSamplerNormalMulti_h