//  matrix_m04.cpp test the matrix channels in c08.cpp.

//  this file has a triple purpose:  1) my testing.  2) demonstration.  3) your computations.



#include <stdio.h>

#include "matrix_dll.h"
#include "StdAfx.h"
//#include "matrix_h0s.h"
//#include "matrix_h01.h"
//#include "matrix_h02.h"
//#include "matrix_h03.h"
//#include "matrix_h04.h"
//#include "matrix_h05.h"
//#include "matrix_h06.h"
//#include "matrix_h07.h"
//#include "matrix_h08.h"
#include "matrix_h09.h"


// get yourself enough arrays.  any unused ones consume neglegible storage.
#define amt_of_arrays 25


matrix_preferences mat_pref;		// declare the global preferences.














int main() {
	FILE *outfile, *datafile, *spherefile, *ellipsoidfile;
	MAT array[amt_of_arrays], scratch[amt_of_arrays];
	CHN chan[amt_of_arrays];
	double residue, residue_identity, norm_average, relative_variance, relative_average, R, Theta, rho;
	double start_LDIAG, ratio_LDIAG, start_MDIAG, ratio_MDIAG, start_RDIAG, ratio_RDIAG, relative_chan[ELEVEN_MX + 1];
	unsigned long i, sw;
	unsigned long counter=0, left_size, right_size, size, multiplier, amount, sort=0, construct;
	unsigned long sw_LBRA, sw_MBRA, sw_MKET, sw_RKET, sw_LAVE, sw_RAVE, answer, method;
	int y=0;
	char title[200];
	extern matrix_preferences mat_pref;		// so that you may access them.

	outfile = fopen("test.txt", "w");

	for(i=0; i<amt_of_arrays; i++) {	// required for the automatic storage allocation to work.
		array[i].mx = 0; scratch[i].mx = 0;
		matrix_initialize_channel(&chan[i]); };



	// *********** the preferences ***************

	matrix_set_preferences();			// required to set the defaults.
//	for safety, do this even if you intend to over-ride all of them.

	// over-ride any (or all) of them.
//	mat_pref.eigenpairs=1;		// 0=system default 1=Jacobi 2=squaring.
	mat_pref.Jacobi_style=0;	// 0=original 1=simplified - by rows  2=simplified - by super-diagonals.
//	mat_pref.test_Jac = 1.0e-18;	// test criterion for Jacobi.
//	mat_pref.Jacobi_angle=0;	// 0=always+ 1=always- 2=small 3=large.  add 4 to have cosine negative.
	mat_pref.random=1;			// 0=system default 1=199017 2=352345383.
	mat_pref.seed=54321;		// the seed for the pseudo-random number generator.

//	mat_pref.test_Jac = 1.0e-10;	// temporarily. <<<<<<<<<<<<<<<<<
//	mat_pref.singular = 1.0e-10;	// temp 


	// ********** ask for the test situation to be specified.  ******************
	
	printf("this is m04\n");
	matrix_about();		// displays the copyright notice.
capo:

	left_size = 5; multiplier = 1000000; amount = 300;
	right_size = left_size; size = left_size + right_size;


	printf("Generating a pseudo-random Gaussian event data-file.\n");
	datafile = fopen("events1.txt", "w");
	y = matrix_write_event_data(datafile, "this_is_a_test", left_size, right_size, multiplier, amount);
	fclose(datafile);
	if(y != 0) goto tie_point;

	printf("Computing its variance & average and assessing the quality of the result.\n");
	datafile = fopen("events1.txt", "r");	//  variance   average
	y = matrix_read_event_data(datafile, title, &array[0], &array[1]);
	fclose(datafile);
	if(y != 0) goto tie_point;
	printf("%s\n", title);
//	matrix_display(&array[0], "variance");
//	matrix_display(&array[1], "average");
	y = matrix_isit_symmetric(&array[0], &residue, &scratch[0], &scratch[1], &scratch[2], &scratch[3]);
	if(y != 0) goto tie_point;
	y = matrix_isit_identity(&array[0], &residue_identity, &scratch[0], &scratch[1], &scratch[2]);
	if(y != 0) goto tie_point;
	y = matrix_norm(&array[1], 1, &norm_average, &scratch[0]);
	if(y != 0) goto tie_point;
	printf("generated sphere:  is the variance symmetric? %e \n", residue);
	printf("   is it spherical? %e  norm of average %e \n\n", residue_identity, norm_average);

	printf("Constructing a pseudo-random Gaussian event data-file.\n"); 
	sort = 0;							//  bra        diag       ket
	y = matrix_square_root(&array[0], sort, &array[2], &array[3], &array[4], &scratch[0],
			&array[5], &scratch[1], &scratch[2]);	// #5 is the correcting array.
	(array[6]).imx = (array[1]).imx; (array[6]).jmx = (array[1]).jmx;
	if(y != 0) goto tie_point;
	y = matrix_zero(&array[6]);
	if(y != 0) goto tie_point;
	matrix_propagate_default(&array[1], &array[6]);
	datafile = fopen("events1.txt", "r");
	spherefile = fopen("eventsO.txt", "w");			//  old average corrting  new average
	y = matrix_read_process_event_data1(datafile, title, &array[1], &array[5], &array[6],
			spherefile, "this_is_a_sphere", &scratch[0]);
	fclose(datafile);
	fclose(spherefile);
	if(y != 0) goto tie_point;

	printf("Computing its variance & average and assessing the quality of the result.\n");
	spherefile = fopen("eventsO.txt", "r");	//    variance   average
	y = matrix_read_event_data(spherefile, title, &array[7], &array[8]);
	fclose(spherefile);
	if(y != 0) goto tie_point;
	printf("%s\n", title);
//	matrix_display(&array[7], "variance");
//	matrix_display(&array[8], "average");
	y = matrix_isit_symmetric(&array[7], &residue, &scratch[0], &scratch[1], &scratch[2], &scratch[3]);
	if(y != 0) goto tie_point;
	y = matrix_isit_identity(&array[7], &residue_identity, &scratch[0], &scratch[1], &scratch[2]);
	if(y != 0) goto tie_point;
	y = matrix_norm(&array[8], 1, &norm_average, &scratch[0]);
	if(y != 0) goto tie_point;
	printf("construced sphere:  is the variance symmetric? %e \n", residue);
	printf("   is it spherical? %e  norm of average %e \n\n", residue_identity, norm_average);

	printf("Constructing a pseudo-random ellipsoidal Gaussian event data-file \n");


//	construction of the variance and the average. =======================================

	method = 0;		// 0=easy way, 1=component way.

	if(method == 0) {
//	the easy way!  ===========================================
		y = matrix_random_positive_symmetric(size, &array[9], &scratch[0]);		// the desired variance.
		if(y != 0) goto tie_point;
		(array[9]).imp = left_size; (array[9]).jmp = left_size; (array[9]).part = 1;
//		matrix_display(&array[9], "the desired variance");
		(array[10]).imx = 1; (array[10]).jmx = size;
		y = matrix_uniformly_random(&array[10]);	// the desired average.
		if(y != 0) goto tie_point;
		(array[10]).jmp = left_size; (array[10]).part = 1;
//		matrix_display(&array[10], "the desired average").
	}

	else {
//	the component way! =======================================

		sw_LBRA = 2;		// 0=identity, 2=random ortho-normal.
		sw_MBRA = 2;		// 0=identity, 2=random ortho-normal.
		sw_MKET = 1;		// 0=identity, 1=transpose of MBRA, 2=random ortho-normal.
		sw_RKET = 1;		// 0=identity, 1=transpose of LBRA, 2=random ortho-normal.
		sw_LAVE = 1;		// 0=zero, 1=random.
		sw_RAVE = 1;		// 0=zero, 1=random.
		start_LDIAG = 1.0;		ratio_LDIAG = 0.6;	// by convention, these should be strictly positive.
		start_MDIAG = 0.9;		ratio_MDIAG = 0.7;	// the start should be in the open interval (-1, 1) and ratio in (-1, 0) union (0, 1).
		start_RDIAG = 1.0;		ratio_RDIAG = 0.5;	// by convention, these should be strictly positive.


		y = matrix_channel_xxx(left_size, sw_LBRA, sw_MBRA, sw_MKET, sw_RKET, sw_LAVE, sw_RAVE,
				start_LDIAG, ratio_LDIAG, start_MDIAG, ratio_MDIAG, start_RDIAG, ratio_RDIAG,
				&chan[1], &scratch[0], &scratch[1], &scratch[2], &scratch[3], &R, &Theta, &rho,
				&scratch[4], &scratch[5], &scratch[6], &scratch[7], &scratch[8], &scratch[9]);
		if(y != 0) goto tie_point;
		matrix_display(&scratch[0], "constructed channel:  exterior bra");
		matrix_display(&scratch[1], "   exterior diag");
		matrix_display(&scratch[2], "   interior bra");
		matrix_display(&scratch[3], "   interior diag");
		printf("   R=%e Theta=%e, rho=%e\n", R, Theta, rho);
		y = matrix_array_bar_channel(&chan[1].eleven[ARR], &chan[1].eleven[AVE], &chan[2],
				&scratch[0], &scratch[1], &scratch[2], &scratch[3], &scratch[4],
				&residue, &scratch[5], &scratch[6], &scratch[7], &scratch[8], &scratch[9],
				&scratch[10], &scratch[11], &scratch[12], &scratch[13], &scratch[14],
				&scratch[15], &scratch[16]);
		if(y != 0) goto tie_point;
		y = matrix_isit_equal_channels(&chan[1], &chan[2], relative_chan, &scratch[0], &scratch[1]);
		if(y != 0) goto tie_point;
		for(i=0; i<ELEVEN_MX; i++) {
			printf("i=%d relative=%e  ", i, relative_chan[i]); };
		printf("\nis the constructed channel valid?  relative analysis %e \n", relative_chan[ELEVEN_MX]);
		y = matrix_copy(&chan[1].eleven[ARR], &array[9]);
		if(y != 0) goto tie_point;
		y = matrix_copy(&chan[1].eleven[AVE], &array[10]);
		if(y != 0) goto tie_point; };


//	the end of the construction of the variance and the average. ========================

//	check the array for positive-definite  ======

//											bra			 diag		  ket
	y = matrix_factor_symmetric(&array[9], &scratch[0], &scratch[1], &scratch[2], 0, &scratch[3]);
	if(y != 0) goto tie_point;
	y = matrix_diagonal_isit_positive(&scratch[1], &answer);
	if(y != 0) goto tie_point;
	printf("is the array positive-definite?  answer=%d \n", answer);

	sort = 0;							//  bra         diag        ket
	y = matrix_square_root(&array[9], sort, &array[11], &array[12], &array[13], &array[14],
			&scratch[0], &scratch[1], &scratch[2]);	// #14 is the driving array.
	spherefile = fopen("eventsO.txt", "r");
	ellipsoidfile = fopen("eventsE.txt", "w");			//  old average driving     new average
	y = matrix_read_process_event_data1(spherefile, title, &array[6], &array[14], &array[10],
			ellipsoidfile, "this_is_an_ellipsoid", &scratch[0]);
	fclose(spherefile);
	fclose(ellipsoidfile);
	if(y != 0) goto tie_point;

	printf("Computing its variance & average and assessing the quality of the result.\n");
	ellipsoidfile = fopen("eventsE.txt", "r");	//    variance   average
	y = matrix_read_event_data(ellipsoidfile, title, &array[15], &array[16]);
	fclose(ellipsoidfile);
	if(y != 0) goto tie_point;
	printf("%s\n", title);
	matrix_display(&array[15], "variance");
	matrix_display(&array[16], "average");
	y = matrix_isit_symmetric(&array[15], &residue, &scratch[0], &scratch[1], &scratch[2], &scratch[3]);
	if(y != 0) goto tie_point;
	y = matrix_isit_equal(&array[9], &array[15], &relative_variance, &scratch[0], &scratch[1]);
	if(y != 0) goto tie_point;
	y = matrix_isit_equal(&array[10], &array[16], &relative_average, &scratch[0], &scratch[1]);
	if(y != 0) goto tie_point;
	printf("constructed ellipsoid:  is the variance symmetric %e \n", residue); 
	printf("   is the variance as desired? %e  is the average as desired? %e \n\n",
			relative_variance, relative_average);

	construct = 1;		// 0=generate, 1=construct.
	printf("Now, let us do it all in one step. construct=%d\n", construct);
	y = matrix_ellipsoidal_events(&array[9], &array[10], construct, multiplier, amount,
			"direct_ellipsoidal.txt", "direct_ellipsoidal",
			&array[15], &array[16], &residue, &relative_variance, &relative_average,
			"scratch_events.txt", "scratch_events",
			&scratch[0], &scratch[1], &scratch[2], &scratch[3], &scratch[4],
			&scratch[5], &scratch[6]);
	if(y != 0) goto tie_point;
	matrix_display(&array[15], "variance");
	matrix_display(&array[16], "average");
	printf("constructed ellipsoid:  is the variance symmetric %e \n", residue); 
	printf("   is the variance as desired? %e  is the average as desired? %e \n\n",
			relative_variance, relative_average); 

	y = matrix_channel_rate(&array[15], &R, &Theta, &rho, &scratch[1], &scratch[5], &scratch[3],
			&scratch[7], &scratch[5], &scratch[6]);
	if(y != 0) goto tie_point;
	printf("R=%e  Theta=%e  rho=%e \n", R, Theta, rho);
	printf("analyzing the channel.\n");
	y = matrix_array_bar_channel(&array[15], &array[16], &chan[0], &scratch[0], &scratch[1],
			&scratch[2], &scratch[3], &scratch[4], &residue, &scratch[5], &scratch[6],
			&scratch[7], &scratch[8], &scratch[9], &scratch[10], &scratch[11], &scratch[12],
			&scratch[13], &scratch[14], &scratch[15], &scratch[16]);
	if(y != 0) goto tie_point;
	printf("residue=%e\n", residue);
	printf("computing the regression.\n");


	for(sw=0; sw<=1; sw++) {
//	sw = 0;		// 0=x to y forward, 1=y to x forward, 2=x to y inverse, 3=y to x inverse.
		y = matrix_events_regress(&chan[0], "direct_ellipsoidal.txt", "regression.txt",
				"this_is_the_regression", "residue.txt", "this_is_the_residue",
				sw, &array[17], &residue, &R, &Theta, &rho,
				 &array[18], &norm_average, &scratch[0], &scratch[1],
				 &scratch[2], &scratch[3], &scratch[4], &scratch[5], &scratch[6]);
		printf("the direct ellipsoidal events file has had its regression computed.\n");
		printf("  symmetric? %e  ave-norm=%e  R=%e Theta=%e rho=%e \n",
			residue, norm_average, R, Theta, rho);
		matrix_display(&array[17], "variance of residue");
		matrix_display(&array[18], "average of residue");
		switch(sw) {
			case 0:
				matrix_display(&chan[0].eleven[VarRESyy], "predicted residue variance yy");
				break;
			case 1:
				matrix_display(&chan[0].eleven[VarRESxx], "predicted residue variance xx");
				break; }; };



	printf("so far, so good. \n");




	counter++;
	goto tie_point;






goto capo;		// del capo al fine.



	// ********* finish.  ************************************8

tie_point:		// fine.
	for(i=0; i<amt_of_arrays; i++) {	// clean-up the storage.
		y += matrix_dealloc(&array[i]) + matrix_dealloc(&scratch[i]);
		y += matrix_channel_dealloc(&chan[i]); };
	fclose(outfile);
	return(0); };





// copyright (c) 2003,4 by R.I. 'Scibor-Marchocki.  last modified Thursday 13-th May 2004.

// eof