Index: /bdm/estim/libPF.h
===================================================================
--- /bdm/estim/libPF.h (revision 33)
+++ /bdm/estim/libPF.h (revision 60)
@@ -19,4 +19,7 @@
 
 using namespace itpp;
+
+//UGLY HACK
+extern double PF_SSAT;//used for StrSim:06 test... if length>0 the value is written.
 
 /*!
@@ -60,5 +63,5 @@
 
 class MPF : public PF {
-	BM_T* Bms[1000];
+	BM_T* Bms[10000];
 
 	//! internal class for MPDF providing composition of eEmp with external components
@@ -105,5 +108,5 @@
 		//
 
-		if ( n>1000 ) it_error ( "increase 1000 here!" );
+		if ( n>10000 ) it_error ( "increase 10000 here!" );
 
 		for ( int i=0;i<n;i++ ) {
@@ -126,4 +129,7 @@
 		for ( int i=0;i<n;i++ ) {Bms[i]->condition ( _samples ( i ) );}
 	}
+
+//SimStr:
+	double SSAT;	
 };
 
@@ -132,10 +138,16 @@
 	int i;
 	vec lls ( n );
+	vec llsP ( n );
 	ivec ind;
 	double mlls=-std::numeric_limits<double>::infinity();
 
+	// StrSim:06
+	double sumLWL=0.0;
+	double sumL2WL=0.0;
+	double WL = 0.0;
+
 	for ( i=0;i<n;i++ ) {
 		//generate new samples from paramater evolution model;
-		_samples ( i ) = par.samplecond ( _samples ( i ), lls ( i ) );
+		_samples ( i ) = par.samplecond ( _samples ( i ), llsP ( i ) );
 		Bms[i]->condition ( _samples ( i ) );
 		Bms[i]->bayes ( dt );
@@ -144,4 +156,13 @@
 		if ( lls ( i ) >mlls ) mlls=lls ( i ); //find maximum likelihood (for numerical stability)
 	}
+	
+	if ( true ) {
+		for ( i=0;i<n;i++ ) {
+			WL = _w ( i ) *exp ( llsP ( i ) ); //using old weights!
+			sumLWL += exp ( lls ( i ) ) *WL;
+			sumL2WL += exp ( 2*lls ( i ) ) *WL;
+		}
+		SSAT  = sumL2WL/(sumLWL*sumLWL);
+	}
 
 	// compute weights
@@ -150,7 +171,9 @@
 	}
 
-	_w /=sum(_w); //?
-
-	double eff = 1.0/( _w*_w );
+	_w /=sum ( _w ); //?
+
+
+
+	double eff = 1.0/ ( _w*_w );
 	if ( eff < ( 0.1*n ) ) {
 		ind = est.resample();
@@ -161,5 +184,5 @@
 				//fixme this would require new assignment operator
 				// *Bms[i] = *Bms[ind ( i ) ];
-				
+
 				// poor-man's solution: replicate constructor here
 				// copied from MPF::MPF
@@ -170,5 +193,5 @@
 			}
 		};
-		cout << '.'; 
+		cout << '.';
 	}
 }
Index: /bdm/estim/libPF.cpp
===================================================================
--- /bdm/estim/libPF.cpp (revision 32)
+++ /bdm/estim/libPF.cpp (revision 60)
@@ -4,5 +4,4 @@
 
 using std::endl;
-
 
 void PF::bayes ( const vec &dt ) {
Index: /bdm/stat/libEF.cpp
===================================================================
--- /bdm/stat/libEF.cpp (revision 32)
+++ /bdm/stat/libEF.cpp (revision 60)
@@ -60,5 +60,5 @@
 
 vec mgamma::samplecond ( vec &cond, double &ll ) {
-	*_beta=k / cond ;
+	this->condition(cond );
 	vec smp = epdf.sample();
 	ll = epdf.evalpdflog ( smp );
Index: /bdm/stat/libEF.h
===================================================================
--- /bdm/stat/libEF.h (revision 51)
+++ /bdm/stat/libEF.h (revision 60)
@@ -96,5 +96,5 @@
 	double eval ( const vec &val ) const ;
 	double evalpdflog ( const vec &val ) const;
-	vec mean()const {return mu;}
+	vec mean() const {return mu;}
 
 //Access methods
@@ -140,5 +140,5 @@
 	//! Returns poiter to alpha and beta. Potentially dengerous: use with care!
 	void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
-	vec mean()const {vec pom(alpha); pom/=beta; return pom;}
+	vec mean() const {vec pom ( alpha ); pom/=beta; return pom;}
 };
 /*
@@ -228,4 +228,5 @@
 */
 class mgamma : public mEF {
+protected:
 	//! Internal epdf that arise by conditioning on \c rvc
 	egamma epdf;
@@ -247,4 +248,31 @@
 };
 
+/*!
+\brief  Gamma random walk around a fixed point
+
+Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, $p$. $k$ is the coefficient of the geometric combimation
+\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
+
+Standard deviation of the random walk is proportional to one $k$-th the mean.
+This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+*/
+class mgamma_fix : public mgamma {
+protected:
+	double l;
+	vec refl;
+public:
+	//! Constructor
+	mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
+	//! Set value of \c k
+	void set_parameters ( double k0 , vec ref0, double l0 ) {
+		mgamma::set_parameters ( k0 );
+		refl=pow ( ref0,1.0-l0 );l=l0;
+	};
+
+	void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
+};
+
 //! Switch between various resampling methods.
 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
@@ -264,5 +292,5 @@
 public:
 	//! Default constructor
-	eEmp ( const RV &rv0 ,int n0) :epdf ( rv0 ),n(n0),w(n),samples(n) {};
+	eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
 	//! Set sample
 	void set_parameters ( const vec &w0, epdf* pdf0 );
@@ -276,7 +304,8 @@
 	vec sample() const {it_error ( "Not implemented" );return 0;}
 	//! inherited operation : NOT implemneted
-	double evalpdflog(const vec &val) const {it_error ( "Not implemented" );return 0.0;}
-	vec mean()const {vec pom=zeros(rv.count()); 
-		for (int i=0;i<n;i++){pom+=samples(i)*w(i);}
+	double evalpdflog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
+	vec mean() const {
+		vec pom=zeros ( rv.count() );
+		for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
 		return pom;
 	}
@@ -287,5 +316,5 @@
 
 template<class sq_T>
-enorm<sq_T>::enorm ( RV &rv ) :eEF(rv), mu ( rv.count() ),R ( rv.count() ),_iR ( rv.count() ),cached ( false ),dim ( rv.count() ) {};
+enorm<sq_T>::enorm ( RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),_iR ( rv.count() ),cached ( false ),dim ( rv.count() ) {};
 
 template<class sq_T>
@@ -320,5 +349,5 @@
 
 template<class sq_T>
-mat enorm<sq_T>::sample ( int N )const {
+mat enorm<sq_T>::sample ( int N ) const {
 	mat X ( dim,N );
 	vec x ( dim );
@@ -346,5 +375,5 @@
 template<class sq_T>
 double enorm<sq_T>::evalpdflog ( const vec &val ) const {
-	if ( !cached ) {it_error("this should not happen, see cached");}
+	if ( !cached ) {it_error ( "this should not happen, see cached" );}
 
 	// 1.83787706640935 = log(2pi)
Index: /tests/pmsm_sim2.cpp
===================================================================
--- /tests/pmsm_sim2.cpp (revision 60)
+++ /tests/pmsm_sim2.cpp (revision 60)
@@ -0,0 +1,203 @@
+/*
+  \file
+  \brief Models for synchronous electric drive using IT++ and BDM
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#include <itpp/itbase.h>
+#include <estim/libKF.h>
+#include <estim/libPF.h>
+#include <stat/libFN.h>
+
+#include "pmsm.h"
+#include "simulator.h"
+
+#include <netcdfcpp.h>
+void write_to_nc ( NcFile &nc, mat &X, std::string Xn, Array<std::string> A ) {
+	char tmpstr[200];
+	int Len = X.rows();
+
+	sprintf ( tmpstr,"%s.length",Xn.c_str() );
+	NcDim* lengt = nc.add_dim ( tmpstr, ( long ) Len );
+	for ( int j=0; j<X.cols(); j++ ) {
+		if ( j<A.length() )
+			sprintf ( tmpstr,"%s_%s",Xn.c_str(), ( A ( j ) ).c_str() );
+		else
+			sprintf ( tmpstr,"%s_%d",Xn.c_str(),j );
+		// Create variables and their attributes
+		NcVar* P = nc.add_var ( tmpstr, ncDouble, lengt );
+		const double* Dp = X._data();
+		P->put ( &Dp[j*Len],Len );
+	}
+}
+
+
+using namespace itpp;
+//!Extended Kalman filter with unknown \c Q
+class EKF_unQ : public EKFCh , public BMcond {
+public:
+	//! Default constructor
+	EKF_unQ ( RV rx, RV ry,RV ru,RV rQ ) :EKFCh ( rx,ry,ru ),BMcond ( rQ ) {};
+	void condition ( const vec &Q0 ) {
+		Q.setD ( Q0,0 );
+		//from EKF
+		preA.set_submatrix ( dimy+dimx,dimy,Q._Ch() );
+	};
+};
+
+void set_simulator_t(double &Ww) {
+
+        if (t>0.0) x[8]=1.2;    // 1A //0.2ZP
+        if (t>0.4) x[8]=10.8;   // 9A
+        if (t>0.6) x[8]=25.2;  // 21A
+
+        if (t>0.7) Ww=2.*M_PI*10.;
+        if (t>1.0) x[8]=1.2;    // 1A
+        if (t>1.2) x[8]=10.8;   // 9A
+        if (t>1.4) x[8]=25.2;  // 21A
+
+        if (t>1.6) Ww=2.*M_PI*50.;
+        if (t>1.9) x[8]=1.2;    // 1A
+        if (t>2.1) x[8]=10.8;   // 9A
+        if (t>2.3) x[8]=25.2;  // 21A
+
+        if (t>2.5) Ww=2.*M_PI*100;
+        if (t>2.8) x[8]=1.2;    // 1A
+        if (t>3.0) x[8]=10.8;   // 9A
+        if (t>3.2) x[8]=25.2;  // 21A
+
+        if (t>3.4) Ww=2.*M_PI*150;
+        if (t>3.7) x[8]=1.2;    // 1A
+        if (t>3.9) x[8]=10.8;   // 9A
+        if (t>4.1) x[8]=25.2;  // 21A
+
+        if (t>4.3) Ww=2.*M_PI*0;
+        if (t>4.8) x[8]=-1.2;    // 1A
+        if (t>5.0) x[8]=-10.8;   // 9A
+        if (t>5.2) x[8]=-25.2;  // 21A
+
+        if (t>5.4) Ww=2.*M_PI*(-10.);
+        if (t>5.7) x[8]=-1.2;    // 1A
+        if (t>5.9) x[8]=-10.8;   // 9A
+        if (t>6.1) x[8]=-25.2;  // 21A
+
+        if (t>6.3) Ww=2.*M_PI*(-50.);
+        if (t>6.7) x[8]=-1.2;    // 1A
+        if (t>6.9) x[8]=-10.8;   // 9A
+        if (t>7.1) x[8]=-25.2;  // 21A
+
+        if (t>7.3) Ww=2.*M_PI*(-100.);
+        if (t>7.7) x[8]=-1.2;    // 1A
+        if (t>7.9) x[8]=-10.8;   // 9A
+        if (t>8.1) x[8]=-25.2;  // 21A
+        if (t>8.3) x[8]=10.8;   // 9A
+        if (t>8.5) x[8]=25.2;  // 21A
+}
+
+int main() {
+	// Kalman filter
+	int Ndat = 30000;
+	double h = 1e-6;
+	int Nsimstep = 125;
+	int Npart = 100;
+	
+	//StrSim:06:
+	vec SSAT(Ndat);
+
+	// internal model
+	IMpmsm fxu;
+	//                  Rs    Ls        dt       Fmag(Ypm) kp   p    J     Bf(Mz)
+	fxu.set_parameters ( 0.28, 0.003465, Nsimstep*h, 0.1989,   1.5 ,4.0, 0.04, 0.0 );
+	// observation model
+	OMpmsm hxu;
+
+	vec mu0= "0.0 0.0 0.0 0.0";
+	vec Qdiag ( "0.01 0.01 0.0001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001
+	vec Rdiag ( "0.005 0.005" ); //var(diff(xth)) = "0.034 0.034"
+	chmat Q ( Qdiag );
+	chmat R ( Rdiag );
+	EKFCh KFE ( rx,ry,ru );
+	KFE.set_parameters ( &fxu,&hxu,Q,R );
+	KFE.set_est ( mu0, chmat ( 1*ones ( 4 ) ) );
+
+	RV rQ ( "100","{Q}","4","0" );
+	EKF_unQ KFEp ( rx,ry,ru,rQ );
+	KFEp.set_parameters ( &fxu,&hxu,Q,R );
+	KFEp.set_est ( mu0, chmat ( 1*ones ( 4 ) ) );
+
+	mgamma_fix evolQ ( rQ,rQ );
+	MPF<EKF_unQ> M ( rx,rQ,evolQ,evolQ,Npart,KFEp );
+	// initialize
+	evolQ.set_parameters ( 1000.0 ,Qdiag, 0.5); //sigma = 1/10 mu
+	evolQ.condition ( Qdiag ); //Zdenek default
+	epdf& pfinit=evolQ._epdf();
+	M.set_est ( pfinit );
+	evolQ.set_parameters ( 100000.0, Qdiag, 0.5 );
+
+	//
+
+	epdf& KFEep = KFE._epdf();
+	epdf& Mep = M._epdf();
+
+	mat Xt=zeros ( Ndat ,9 ); //true state from simulator
+	mat Dt=zeros ( Ndat,4+2 ); //observation
+	mat XtE=zeros ( Ndat, 4 );
+	mat XtM=zeros ( Ndat,4+4 ); //Q + x
+
+	// SET SIMULATOR
+	pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
+	double Ww=0.0;
+	vec dt ( 2 );
+	vec ut ( 2 );
+
+	for ( int tK=1;tK<Ndat;tK++ ) {
+		//Number of steps of a simulator for one step of Kalman
+		for ( int ii=0; ii<Nsimstep;ii++ ) {
+			//simulator
+			set_simulator_t(Ww);
+			pmsmsim_step ( Ww );
+		};
+		// collect data
+		ut ( 0 ) = KalmanObs[0];
+		ut ( 1 ) = KalmanObs[1];
+		dt ( 0 ) = KalmanObs[2];
+		dt ( 1 ) = KalmanObs[3];
+
+		//estimator
+		KFE.bayes ( concat ( dt,ut ) );
+		M.bayes ( concat ( dt,ut ) );
+		SSAT(tK) = M.SSAT;
+		
+		Xt.set_row ( tK,vec ( x,9 ) ); //vec from C-array
+		Dt.set_row ( tK, concat ( dt,ut,vec_1(sqrt(pow(ut(0),2)+pow(ut(1),2))), vec_1(sqrt(pow(dt(0),2)+pow(dt(1),2))) ) );
+		XtE.set_row ( tK,KFEep.mean() );
+		XtM.set_row ( tK,Mep.mean() );
+	}
+
+	it_file fou ( "pmsm_sim.it" );
+
+	fou << Name ( "xth" ) << Xt;
+	fou << Name ( "Dt" ) << Dt;
+	fou << Name ( "xthE" ) << XtE;
+	fou << Name ( "xthM" ) << XtM;
+	fou << Name ( "SSAT" ) << SSAT;
+	//Exit program:
+
+	////////////////
+	// Just Testing
+	NcFile nc ( "pmsm_sim.nc", NcFile::Replace ); // Create and leave in define mode
+	if ( ! nc.is_valid() ) {	std::cerr << "creation of NCFile failed."<<endl;}
+
+	write_to_nc ( nc,Xt,"X","{isa isb om th }" );
+	write_to_nc ( nc,XtM,"XtM","{q1 q2 isa isb om th }" );
+	write_to_nc ( nc,XtE,"XE","{isa isb om th }" );
+	write_to_nc ( nc,Dt,"Dt","{isa isb ua ub }" );
+	return 0;
+}
+
Index: /tests/CMakeLists.txt
===================================================================
--- /tests/CMakeLists.txt (revision 54)
+++ /tests/CMakeLists.txt (revision 60)
@@ -46,3 +46,5 @@
 link_directories (${BDM_BINARY_DIR}/simulator_zdenek)
 add_executable (pmsm_sim pmsm_sim.cpp)
+add_executable (pmsm_sim2 pmsm_sim2.cpp)
 target_link_libraries (pmsm_sim ${BdmLibs} pmsmsim netcdf_c++)
+target_link_libraries (pmsm_sim2 ${BdmLibs} pmsmsim netcdf_c++)
Index: /BDM.kdevelop.filelist
===================================================================
--- /BDM.kdevelop.filelist (revision 45)
+++ /BDM.kdevelop.filelist (revision 60)
@@ -28,4 +28,5 @@
 tests/CMakeLists.txt
 tests/pmsm_sim.cpp
+tests/pmsm_sim2.cpp
 tests/pmsm_unkQ.cpp
 tests/pmsm_unkQpf.cpp
@@ -38,2 +39,3 @@
 tests/testResample.cpp
 tests/testSmp.cpp
+tests/testbidiff.cpp
Index: /BDM.kdevelop
===================================================================
--- /BDM.kdevelop (revision 45)
+++ /BDM.kdevelop (revision 60)
@@ -24,5 +24,5 @@
     <filelistdirectory>/home/smidl/work/mixpp</filelistdirectory>
     <run>
-      <mainprogram>/home/smidl/work/mixpp/tests/pmsm_sim</mainprogram>
+      <mainprogram>/home/smidl/work/mixpp/tests/pmsm_sim2</mainprogram>
       <directoryradio>custom</directoryradio>
       <customdirectory>/home/smidl/work/mixpp</customdirectory>
@@ -171,5 +171,5 @@
     <tree>
       <hidepatterns>*.o,*.lo,CVS,*~,cmake*,Makefile</hidepatterns>
-      <hidenonprojectfiles>true</hidenonprojectfiles>
+      <hidenonprojectfiles>false</hidenonprojectfiles>
       <showvcsfields>false</showvcsfields>
     </tree>
