// Port of Bock-Bargmann R to C++
// dt 7/2004 CodeWarrior 8.3 / OS X

#include <iostream>
#include <ctime> 		 // for timing routines
#include "BockBargmann.h"

int main() {
	using namespace std;
	clock_t time1, time2;
	
	int nvar = 6;
	int npar = 7;
	double N2 = 152/2;
	SymmetricMatrix S(nvar);
	Matrix K(nvar,nvar);
	SymmetricMatrix Hinv(npar);
	
	Real s[] = { 521,
                 477, 576,
                 484, 536, 601,
                 510, 575, 593, 755,
                 523, 580, 598, 718, 797, 
                 528, 584, 613, 722, 751, 802 };
    S << s;
    cout << "Observed covariance matrix:" << endl;
    cout << setprecision(0) << S;
    
	Real k[] = {1,0,0,0,0,0,
				1,1,0,0,0,0,
				1,1,1,0,0,0,
				1,1,1,1,0,0,
				1,1,1,1,1,0,
				1,1,1,1,1,1};
	K << k;
    cout << "K matrix:" << endl;
    cout << setprecision(0) << K;
	cout << setprecision(6);

	ColumnVector p(npar);
	ColumnVector para(npar);
	Real startvalues[] = {500, 60, 20, 80, 20, 25, 50};
	p << startvalues;
	
	BB_LL BB(S, K, p, N2, nvar, npar); // loglikehood "object"

	MLE_D_FI mle_d_fi(BB,20,0.0000001);    // mle "object" 

	para = p;
	time1 = clock();
	mle_d_fi.Fit(para);                       // do the fit
	time2 = clock();
	cout << setprecision(4) << resetiosflags(ios::right);
	cout << "Time required: " << float(time2-time1)/CLOCKS_PER_SEC << " seconds" << endl;
	cout << "Time units: " << CLOCKS_PER_SEC << " ticks/second" << endl;

	cout << setprecision(3) << endl;
	cout << "Estimated parameters:" << endl << para << endl;
	Hinv = BB.FI().i();
	cout << setprecision(1);
	cout << "Inverse Hessian:" << endl << Hinv << endl;
	cout << "Standard Errors:" << setprecision(5) << endl;
	for (int i=1; i <= npar; i++) cout << sqrt(Hinv(i,i)) << endl;

	return 0;
}