/*!
  \file
  \brief TR 2525 file for testing Toy Problem of mpf for Covariance Estimation
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/



#include <estim/libPF.h>
#include <estim/ekf_templ.h>
#include <stat/libFN.h>

#include <stat/loggers_ui.h>
#include <stat/libEF_ui.h>

#include "../pmsm.h"
#include "simulator.h"
#include "../sim_profiles.h"

using namespace bdm;

int main ( int argc, char* argv[] ) {
	const char *fname;
	if ( argc>1 ) {fname = argv[1];	}
	else { fname = "unitsteps.cfg"; }
	UIFile F ( fname );

	int Ndat;
	int Npart;
	double h = 1e-6;
	int Nsimstep = 125;

	vec Qdiag;
	vec Rdiag;

//	mpdf* evolQ ;
	try {
		// Kalman filter
		F.lookupValue ( "ndat", Ndat );
		F.lookupValue ( "Npart",Npart );

//		UIbuild ( F.lookup ( "Qrw" ),evolQ );
		Qdiag= getvec ( F.lookup ( "dQ" ) ); //( "1e-6 1e-6 0.001 0.0001" ); //zdenek: 0.01 0.01 0.0001 0.0001
		Rdiag=getvec ( F.lookup ( "dR" ) );// ( "1e-8 1e-8" ); //var(diff(xth)) = "0.034 0.034"
	}
	catch UICATCH;
// 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";
	chmat Q ( Qdiag );
	chmat R ( Rdiag );
	EKFCh KFE ;
	KFE.set_parameters ( &fxu,&hxu,Q,R );
	KFE.set_est ( mu0, chmat ( zeros ( 4 ) ) );
	KFE.set_rv ( rx );

	RV rQ ( "{Q }","16" );
	EKFCh_chQ KFEp ;
	KFEp.set_parameters ( &fxu,&hxu,Q,R );
	KFEp.set_est ( mu0, chmat ( zeros ( 4 ) ) );

	rwiWishartCh* evolQw = new rwiWishartCh; 
	evolQw->set_parameters(4, 0.1, sqrt(Qdiag),0.99);
	MPF<EKFCh_chQ> M;
	M.set_parameters ( evolQw,evolQw,Npart );
	// initialize
	chmat Ch0(diag(Qdiag));
	evolQw->condition ( vec(Ch0._Ch()._data(),16) ); //Zdenek default
	M.set_statistics ( evolQw->_e() , &KFEp );
	//

	M.set_rv ( concat ( rQ,rx ) );

	dirfilelog *L; UIbuild ( F.lookup ( "logger" ), L );// ( "exp/mpf_test",100 );
	int l_X = L->add ( rx, "xt" );
	int l_D = L->add ( concat ( ry,ru ), "" );
	int l_Q= L->add ( rQ, "" );
	int l_fullQ= L->add ( rQ, "full" );

	KFE.set_options ( "logbounds" );
	KFE.log_add ( L,"KF" );
	M.set_options ( "logbounds" );
	M.log_add ( L,"M" );
	L->init();

	// SET SIMULATOR
	pmsmsim_set_parameters ( 0.28,0.003465,0.1989,0.0,4,1.5,0.04, 200., 3e-6, h );
	vec dt ( 2 );
	vec ut ( 2 );
	vec xt ( 4 );
	vec xtm=zeros ( 4 );
	double Ww=0.0;
	vec vecW=getvec ( F.lookup ( "profile" ) );
	
	mat tQ=diag(Qdiag);
	mat tChQ=chol(tQ);

	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
			sim_profile_vec01t ( Ww,vecW );
			pmsmsim_step ( Ww );
		};
		ut ( 0 ) = KalmanObs[4];
		ut ( 1 ) = KalmanObs[5];
		xt = fxu.eval ( xtm,ut ) + tChQ.T() *randn ( 4 );
		dt = hxu.eval ( xt,ut );
		xtm = xt;

		//Variances
/*		if ( tK==1000 )  tQ ( 0,0 ) *=10; 
		if ( tK==2000 ) tQ ( 0,0 ) /=10;
		if ( tK==3000 )  tQ( 1,1 ) *=10;
		if ( tK==4000 ) tQ( 1,1 ) /=10;
		if ( tK==5000 )  tQ( 2,2 ) *=10;
		if ( tK==6000 ) tQ( 2,2 ) /=10;
		if ( tK==7000 )  tQ( 3,3 ) *=10;
		if ( tK==8000 ) tQ( 3,3 ) /=10;*/
		
		if (tK>1000) {tQ(0,1)=0.5*sqrt(tQ(0,0)*tQ(1,1));tQ(1,0)=tQ(0,1);}
		if (tK>2000) {tQ(0,1)=0; tQ(1,0)=tQ(0,1);}
		
		if (tK>3000) {tQ(2,3)=-0.5*sqrt(tQ(2,2)*tQ(3,3)); tQ(3,2)=tQ(2,3);}
		if (tK>4000) {tQ(2,3)=0; tQ(3,2)=tQ(2,3);}
		
		if (tK>5000) {tQ(0,2)=0.9*sqrt(tQ(0,0)*tQ(2,2)); tQ(2,0)=tQ(0,2);}
		if (tK>6000) {tQ(0,2)=0; tQ(2,0)=tQ(0,2);}

		tChQ=chol(tQ);
		
		//estimator
		KFE.bayes ( concat ( dt,ut ) );
		M.bayes ( concat ( dt,ut ) );

		L->logit ( l_X,xt );
		L->logit ( l_D,concat ( dt,ut ) );
		mat Q=diag(Qdiag);
		L->logit ( l_Q,vec(tQ._data(),16) );
		
		mat chQ(4,4);
		copy_vector(16,M._e()->mean()._data(),chQ._data());
		mat fQ=chQ.T()*chQ;
		L->logit ( l_fullQ,vec(fQ._data(),16) );

		KFE.logit ( L );
		M.logit ( L );
		L->step();
	}
	L->finalize();
	//Exit program:

	delete L;
	return 0;
}