/*
	esymrand.c
	Generate random orientations or randomly choose one
	of the symmetry-related orientations.
	Authors: Eduardo Sanz-Garcia and David Belnap
	Created: 20051011          Modified: 20090827
*/

#include "utilities.h"
#include "stdio.h"
#include "options.h"
#include "rwmg.h"
#include "linked_list.h"
#include "symmetry.h"
#include "random_numbers.h"


// Declaration of global variables
extern int 	verbose;		// Level of output to the screen


// Function prototype
int			project_random_views(Bproject* project);
int			project_random_symmetry_related_views(Bsymmetry* sym, Bproject* project);
int			project_change_views_to_asymmetric_unit(Bsymmetry* sym, Bproject* project);


// Usage assistance
const char* use[] = {
" ",
"Usage: esymrand [options] input.star [input.star]",
"----------------------------------------------------",
"Generate random orientations or randomly choose one",
"of the symmetry-related orientations.",
" ",
"Actions:",
"-random                  Produces random orientations inside in the C1 asymmetric unit.",
"-symrand                 Generates all symmetry-related orientations applying the symmetry of",
"                         a point-group and randomly choose one of the symmetry-related views.",
" ",
"Options:",
"-verbose 7               Verbosity of output",
"-symmetry C5             Point group symmetry (default C1).",
"-asu                     Change to orientations to fit into the asymmetric unit.",
" ",
"Output:",
"-output file.star        Output parameter file name.",
" ",
NULL
};


int		main(int argc, char **argv)
{

	Bproject*		project=NULL;		// Data for all micrographs
	Bsymmetry*		sym;				// Symmetry structure	
	char*			outfile = NULL;		// Output STAR file
	Bstring			symmetry_string;	// No symmetry specified
	int			 	random = 0;			// Flag to generate random views
	int				symrand = 0;		// Flag to generate symmetry-related-random views
	int				asu = 0;			// Flag to change views to ASU
	int				i,p;				// i counts the numbers of symmetry-related views
	int				optind;				// Option index
	Option*			option;
	Option*			curropt;

	option = get_option_list(use, argc, argv, &optind);
	for ( curropt = option; curropt; curropt = curropt->next ) {
		if ( strcmp(curropt->tag, "random") == 0 )
			random = 1;
		if ( strcmp(curropt->tag, "symrand") == 0)
			symrand = 1;
		if ( strcmp(curropt->tag, "symmetry") == 0 )
			symmetry_string = get_option_symmetry(curropt->value);
		if ( strcmp(curropt->tag, "asu") == 0 )
		 	asu = 1;
		if ( strcmp(curropt->tag, "output") == 0 )
				outfile = get_option_filename(curropt->value);
	}
	option_kill(option);

	// Errors
	if (random &&  symrand){
		fprintf(stderr, "Error:  Incompatible options!\n\n");
		bexit(-1);
	}

	// Read all the parameter files
	Bstring*			file_list = NULL;
	while ( optind < argc ) string_add(&file_list, argv[optind++]);
	if ( !file_list ) {
		fprintf(stderr, "Error: No parameter files specified!\n");
		bexit(-1);
	}

	// Errors reading STAR file
	project = read_project(file_list);
	if ( project == NULL )  {
		fprintf(stderr, "Error: No input file read!\n");
		bexit(-1);
	}

	// Set sym structure
	sym = init_point_group_symmetry(symmetry_string);

	// Procede with the STAR file
	if (random)
		project_random_views(project);
	if(symrand)
		project_random_symmetry_related_views(sym, project);
	if (asu)
		project_change_views_to_asymmetric_unit(sym, project);

	// Output parameters
	if ( project && outfile ) {
		project_update_comment(project, argc, argv);
		write_project(outfile, project);
	}

	// Memory cleanup
	kill_symmetry(sym);
	project_kill(project);
	bexit(0);

}   // End of main program

/************************************************************************
@Function: project_random_views
@Description:
	Subtitute the original orientations by random views.
@Algorithm:
	Same algorithm than random_views (util/views.c).
@Arguments:
	Bproject* project	project parameter structure.
@Returns:
	int			0.
*************************************************************************/
int			project_random_views(Bproject* project)
{
	if ( !project ) return(0);
	
	Bfield* 			field;
	Bmicrograph*		mg;
	Breconstruction*	rec;
	Bparticle*			part;
	double				irm = 1.0/get_rand_max(), irm2 = 2*irm;
	View*				v = NULL;
	v = (View *) add_item((char **) &v, sizeof(View));

	srandom(getpid());	 // Seed random generator = process ID
	
	if ( verbose & VERB_PROCESS ) {
		printf("Introducing random views:\n");
	}
	
	if ( project->select ) {
		for ( rec = project->rec; rec; rec = rec->next ) {
			for ( part = rec->part; part; part = part->next ) {
				if (verbose){
					printf ("Original view vector: ");
					show_view(part->view);
				}
				v->x = random()*irm2 - 1;
				v->y = random()*irm2 - 1;
				v->z = random()*irm2 - 1;
				v->a = M_PI*(random()*irm2 - 1);
				view_normalize(v);
				part->view=*v;
				if (verbose){
					printf ("Random view vector: ");
					show_view(part->view);
				}
			}
		}
	} else {
		for ( field = project->field; field; field = field->next ) {
			for ( mg = field->mg; mg; mg = mg->next ) {
				for ( part = mg->part; part; part = part->next ) {
					if (verbose){
						printf ("Original view vector: ");
						show_view(part->view);
					}
					v->x = random()*irm2 - 1;
					v->y = random()*irm2 - 1;
					v->z = random()*irm2 - 1;
					v->a = M_PI*(random()*irm2 - 1);
					view_normalize(v);
					part->view=*v;
					if (verbose){
						printf ("Random view vector: ");
						show_view(part->view);
					}
				}
			}
		}
	}

	return(0);
}

/************************************************************************
@Function: project_random_symmetry_related_views
@Description:
	Subtitute the original orientations by one of the symmetry-related views.
@Algorithm:
@Arguments:
	Bsymmetry*	 sym	symmetry structure.
	Bproject* project	project parameter structure.
@Returns:
	int			0.
*************************************************************************/
int			project_random_symmetry_related_views(Bsymmetry* sym, Bproject* project)
{
	if ( !project ) return(0);
	
	Bfield* 			field;
	Bmicrograph*		mg;
	Breconstruction*	rec;
	Bparticle*			part;
	View*				v;
	int			 		nviews = 0;	// Number of symmetry-related views
	double				r;			// Stores a random value in range [0,1) or [0,n)
	int					y;			// Stores a random integer in range [0,n-1]
	int					i;

	srandom(getpid());	// Seed random generator = process ID
	
	if ( verbose & VERB_PROCESS ) {
		printf("Substituting original orientation by randomly assigning one of the symmetry-related views:\n");
		printf(" Symmetry: %s\n", sym->label);
	}
	
	if ( project->select ) {
		for ( rec = project->rec; rec; rec = rec->next ) {
			for ( part = rec->part; part; part = part->next ) {
				if (verbose){
					printf ("Original view vector: ");
					show_view(part->view);
				 }

				v = symmetry_get_all_views(sym, part->view); // Get symmetry-related views
				if (verbose){
					printf ("Symmetry-related views:\n");
					show_views(v);
				}
				
				 // Count the numbers of symmetry-related views		
				 if (!nviews){
					View*	view=v;
					for (nviews=0; view != NULL ; nviews++, view=view->next)
						;
				 }

				r = (double) random() /((double) RAND_MAX + (double) 1); // r is a random floating between [0,1)
				r = r * nviews;	// r is now a random floating point value in the range [0,n)
				y = (int) r;	// y is a random integer in the range [0,n-1]

				for ( i=0; i<y ; i++, v=v->next)
					;
				part->view.x=v->x;
				part->view.y=v->y;
				part->view.z=v->z;
				part->view.a=v->a;
				
				if (verbose){
					printf ("View chosen #%i: ", i);
					show_view(part->view);
				}
			}
		}
	} else {
		for ( field = project->field; field; field = field->next ) {
			for ( mg = field->mg; mg; mg = mg->next ) {
				for ( part = mg->part; part; part = part->next ) {
					if (verbose){
						printf ("Original view vector: ");
						show_view(part->view);
					 }

					v = symmetry_get_all_views(sym, part->view); // Get symmetry-related views
					if (verbose){
						printf ("Symmetry-related views:\n");
						show_views(v);
					}
				
					 // Count the numbers of symmetry-related views		
					 if (!nviews){
						View*	view=v;
						for (nviews=0; view != NULL ; nviews++, view=view->next)
							;
					 }

					r = (double) random() /((double) RAND_MAX + (double) 1); // r is a random floating between [0,1)
					r = r * nviews;	// r is now a random floating point value in the range [0,n)
					y = (int) r;	// y is a random integer in the range [0,n-1]

					for ( i=0; i<y ; i++, v=v->next)
						;
					part->view.x=v->x;
					part->view.y=v->y;
					part->view.z=v->z;
					part->view.a=v->a;
					
					if (verbose){
						printf ("View chosen #%i: ", i);
						show_view(part->view);
					}
				}
			}
		}
	}

	return(0);
}

/************************************************************************
@Function: project_change_views_to_asymmetric_unit
@Description:
	Subtitute the original orientations by one view in the asymmetric unit.
@Algorithm:
@Arguments:
	Bsymmetry*	 sym	symmetry structure.
	Bproject* project	project parameter structure.
@Returns:
	int			0.
*************************************************************************/
int			project_change_views_to_asymmetric_unit(Bsymmetry* sym, Bproject* project)
{
	if ( !project ) return(0);
	
	Bfield* 			field;
	Bmicrograph*		mg;
	Breconstruction*	rec;
	Bparticle*			part;
	
	if ( verbose & VERB_PROCESS ) {
		printf("Shifting original orientation to the asymmetric unit: %s\n", sym->label);
	}
	
	if ( project->select ) {
		for ( rec = project->rec; rec; rec = rec->next ) {
			for ( part = rec->part; part; part = part->next ) {
				change_views_to_asymmetric_unit(sym, &part->view);
				if (verbose){
					printf ("View shifted to ASU: ");
					show_view(part->view);
				}
			}
		}
	} else {
		for ( field = project->field; field; field = field->next ) {
			for ( mg = field->mg; mg; mg = mg->next ) {
				for ( part = mg->part; part; part = part->next ) {
					change_views_to_asymmetric_unit(sym, &part->view);
					if (verbose){
						printf ("View shifted to ASU: ");
						show_view(part->view);
					}
				}
			}
		}
	}

	return(0);
}