Index: library/bdm/itpp_ext.h
===================================================================
--- library/bdm/itpp_ext.h (revision 328)
+++ library/bdm/itpp_ext.h (revision 404)
@@ -26,6 +26,10 @@
 	bvec operator& ( const bvec &a, const bvec &b );
 	bvec operator| ( const bvec &a, const bvec &b );
+	
+	bvec operator>(const vec &t1, const vec &t2); 
 
-// template<class Num_T>
+	bvec operator<(const vec &t1, const vec &t2); 
+
+	// template<class Num_T>
 // void set_subvector(vec &ov, ivec &iv, const Vec<Num_T> &v);
 
@@ -38,4 +42,7 @@
 	}
 
+	const double inf =std::numeric_limits<double>::infinity();
+	
+//#if 0
 	/*!
 	  \brief Gamma distribution
@@ -74,5 +81,5 @@
 	};
 	bool qr ( const mat &A, mat &R );
-
+//#endif
 	//! reimplementation of matlab num2str
 	std::string num2str(double d);
Index: library/bdm/base/bdmbase.h
===================================================================
--- library/bdm/base/bdmbase.h (revision 395)
+++ library/bdm/base/bdmbase.h (revision 404)
@@ -265,4 +265,5 @@
   virtual mat sample_m(int N) const;
   //! Compute log-probability of argument \c val
+  //! In case the argument is out of suport return -Infinity
   virtual double evallog(const vec &val) const {it_error("not implemneted"); return 0.0;};
   //! Compute log-probability of multiple values argument \c val
@@ -381,5 +382,5 @@
     this->condition(cond);
     tmp = ep->evallog(dt);
-    it_assert_debug(std::isfinite(tmp), "Infinite value");
+ //   it_assert_debug(std::isfinite(tmp), "Infinite value");
     return tmp;
   };
Index: library/bdm/stat/merger.cpp
===================================================================
--- library/bdm/stat/merger.cpp (revision 399)
+++ library/bdm/stat/merger.cpp (revision 404)
@@ -5,39 +5,54 @@
 namespace bdm
 {
-	vec merger_base::merge_points ( mat &lW )
-	{
+	vec merger_base::merge_points ( mat &lW ) {
 		int nu=lW.rows();
-		switch ( METHOD )
-		{
+		vec result;
+		ivec indW;
+		bool infexist;
+		switch ( METHOD ) {
 			case ARITHMETIC:
-				return log ( sum ( exp ( lW ) ) ); //ugly!
+				result= log ( sum ( exp ( lW ) ) ); //ugly!
 				break;
 			case GEOMETRIC:
-				return sum ( lW ) /nu;
+				result= sum ( lW ) /nu;
 				break;
 			case LOGNORMAL:
-				vec mu = sum ( lW ) /nu; //mean of logs
-				// vec lam = sum ( pow ( lW,2 ) )-nu*pow ( mu,2 ); ======= numerically unsafe!
-				vec lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
+				vec sumlW=sum ( lW ) ;
+				indW=find((sumlW<inf) & (sumlW>-inf));
+				infexist=(indW.size()<lW.cols());
+				vec mu;
+				vec lam;
+				if (infexist){
+					mu = sumlW(indW) /nu; //mean of logs
+					// 
+					mat validlW=lW.get_cols(indW);
+					lam = sum ( pow ( validlW-outer_product ( ones ( validlW.rows() ),mu ),2 ) );
+				}
+				else {
+					mu = sum ( lW ) /nu; //mean of logs
+					lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
+				}
+				//
 				double coef=0.0;
 				vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere
-				switch ( nu )
-				{
+				switch ( nu ) {
 					case 2:
 						coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) );
-						return  coef*sq2bl + mu ;
+						result =coef*sq2bl + mu ;
 						break;
-					case 3://Ratio of Bessel
-						coef = sqrt ( ( 3*beta-2 ) /3*beta );
-						return log ( besselk ( 0,sq2bl*coef ) ) - log ( besselk ( 0,sq2bl ) ) +  mu;
-						break;
-					case 4:
-						break;
-					default: // Approximate conditional density
+					// case 4: == can be done similar to case 2 - is it worth it???
+					default: // see accompanying document merge_lognorm_derivation.lyx
+						coef = sqrt(1-(nu+1)/(2*beta*nu));
+						result =log(besselk((nu-3)/2, sq2bl*coef))-log(besselk((nu-3)/2, sq2bl)) + mu;
 						break;
 				}
 				break;
 		}
-		return vec ( 0 );
+		if (infexist){
+			vec tmp =-inf*ones(lW.cols());
+			set_subvector(tmp, indW, result);
+			return tmp;
+		}
+		else {return result;}
 	}
 
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 395)
+++ library/bdm/stat/exp_family.h (revision 404)
@@ -51,5 +51,5 @@
 				double tmp;
 				tmp= evallog_nn ( val )-lognc();
-				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 
+// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 
 				return tmp;}
 			//!Evaluate normalized log-probability for many samples
@@ -280,5 +280,6 @@
 			double evallog_nn ( const vec &val ) const
 			{
-				double tmp; tmp= ( beta-1 ) *log ( val );		it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+				double tmp; tmp= ( beta-1 ) *log ( val );
+// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
 				return tmp;
 			};
@@ -290,5 +291,5 @@
 				for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
 				tmp= lgb-lgamma ( gam );
-				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
 				return tmp;
 			};
@@ -489,5 +490,8 @@
 
 			double eval ( const vec &val ) const  {return nk;}
-			double evallog ( const vec &val ) const  {return lnk;}
+			double evallog ( const vec &val ) const  {
+				if (any(val<low) && any(val>high)) {return inf;}
+				else return lnk;
+			}
 			vec sample() const
 			{
Index: library/bdm/stat/merger.h
===================================================================
--- library/bdm/stat/merger.h (revision 399)
+++ library/bdm/stat/merger.h (revision 404)
@@ -242,5 +242,5 @@
 			vec mea = mean();
 
-			cout << sum (w) << "," << w*w << endl;
+// 			cout << sum (w) << "," << w*w << endl;
 
 			mat Tmp = zeros (dim, dim);
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 395)
+++ library/bdm/stat/exp_family.cpp (revision 404)
@@ -24,5 +24,5 @@
 	if ( dimx==1 ) { //same as the following, just quicker.
 		double r = val ( vend ); //last entry!
-		if (r<0) return -1e+100;
+		if (r<0) return -inf;
 		vec Psi ( nPsi+dimx );
 		Psi ( 0 ) = -1.0;
@@ -36,5 +36,5 @@
 		fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) );
 		double ldetR=R.logdet();
-		if (ldetR) return -1e+100;
+		if (ldetR) return -inf;
 		mat Tmp=concat_vertical ( -eye ( dimx ),Th );
 		fsqmat iR ( dimx );
@@ -190,6 +190,6 @@
 	int i;
 
-	if (any(val<=0.)) return -1e100;
-	if (any(beta<=0.)) return -1e100;
+	if (any(val<=0.)) return -inf;
+	if (any(beta<=0.)) return -inf;
 	for ( i=0; i<dim; i++ ) {
 		res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );
Index: library/bdm/itpp_ext.cpp
===================================================================
--- library/bdm/itpp_ext.cpp (revision 343)
+++ library/bdm/itpp_ext.cpp (revision 404)
@@ -63,29 +63,4 @@
 
 // Gamma
-
-	Gamma_RNG::Gamma_RNG ( double a, double b ) {
-		setup ( a,b );
-	}
-
-
-	bvec operator& ( const bvec &a, const bvec &b ) {
-		it_assert_debug ( b.size() ==a.size(), "operator&(): Vectors of different lengths" );
-
-		bvec temp ( a.size() );
-		for ( int i = 0;i < a.size();i++ ) {
-			temp ( i ) = a ( i ) & b ( i );
-		}
-		return temp;
-	}
-
-	bvec operator| ( const bvec &a, const bvec &b ) {
-		it_assert_debug ( b.size() !=a.size(), "operator&(): Vectors of different lengths" );
-
-		bvec temp ( a.size() );
-		for ( int i = 0;i < a.size();i++ ) {
-			temp ( i ) = a ( i ) | b ( i );
-		}
-		return temp;
-	}
 #define log std::log
 #define exp std::exp
@@ -93,4 +68,45 @@
 #define R_FINITE std::isfinite
 
+	bvec operator>(const vec &t1, const vec &t2) {
+		it_assert_debug(t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors");
+		bvec temp(t1.length());
+		for (int i = 0; i < t1.length(); i++)
+			temp(i) = (t1[i] > t2[i]);
+		return temp;
+	}
+	
+	bvec operator<(const vec &t1, const vec &t2) {
+		it_assert_debug(t1.length() == t2.length(), "Vec<>::operator>(): different size of vectors");
+		bvec temp(t1.length());
+		for (int i = 0; i < t1.length(); i++)
+			temp(i) = (t1[i] < t2[i]);
+		return temp;
+	}
+
+
+	bvec operator& ( const bvec &a, const bvec &b ) {
+		it_assert_debug ( b.size() ==a.size(), "operator&(): Vectors of different lengths" );
+
+		bvec temp ( a.size() );
+		for ( int i = 0;i < a.size();i++ ) {
+			temp ( i ) = a ( i ) & b ( i );
+		}
+		return temp;
+	}
+
+	bvec operator| ( const bvec &a, const bvec &b ) {
+		it_assert_debug ( b.size() !=a.size(), "operator&(): Vectors of different lengths" );
+
+		bvec temp ( a.size() );
+		for ( int i = 0;i < a.size();i++ ) {
+			temp ( i ) = a ( i ) | b ( i );
+		}
+		return temp;
+	}
+
+//#if 0
+	Gamma_RNG::Gamma_RNG ( double a, double b ) {
+	setup ( a,b );
+}
 	double  Gamma_RNG::sample() {
 		//A copy of rgamma code from the R package!!
@@ -290,5 +306,5 @@
 	}
 
-
+//#endif
 	std::string num2str(double d){
 		char tmp[20];//that should do
