/* Copyright (C) 2003-2006 Damir Zucic */

/*=============================================================================

				group.c

Purpose:
	Execute group command:  catch a set of macromolecular complexes,
	releasing some others. A single rotation center will be prepared
	for the entire group,  except if grouping was done using keyword
	ALL.  In that case, each macromolecular center will have its own
	rotation center  (the old one).  All atoms,  including the atoms
	which are not selected, will be taken into account. To break the
	group, use the command CATCH.  If keyword ALL is not used,  then
	two or more  complex identifiers  should be specified as command
	arguments. The keyword SAVE may be used to save the whole group,
	where  each structure  will be treated  as  a single chain.  The
	keyword SSC may be used to save the structure as a single chain,
	with  no chain  identifier.  If any structure consist  of two or
	more chains,  the original chain identifiers  will be lost.  The
	keyword  ALL  may be used to group all loaded structures,  where
	each of these structures will retain its own rotation center.

Input:
	(1) Pointer to MolComplexS structure,  with macromol. complexes.
	(2) Number of macromolecular complexes.
	(3) Pointer to RuntimeS structure, with some runtime data.
	(4) Pointer to ConfigS structure, with configuration data.
	(5) Pointer to GUIS structure, with GUI data.
	(6) Pointer to NearestAtomS structure.
	(7) The number of pixels in the main window free area.
	(8) Pointer to refreshI.
	(9) The string which contains  one or more the macromol. complex
	    identifiers.

Output:
	(1) The catch flag (catchF) set for each macromolecular complex.
	(2) The edit mode index will be reset to zero.
	(3) The group flag in RuntimeS structure will be set to one.
	(4) Return value.

Return value:
	(1) Positive (command) code on success.
	(2) Negative (error) code on failure.

Notes:
	(1) This command takes one of more arguments.  To catch only one
	    macromolecular complex, use the command catch.

	(2) This function resets the edit mode index.

========includes:============================================================*/

#include <stdio.h>

#include <string.h>
#include <stdlib.h>

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <X11/Xatom.h>

#include "defines.h"
#include "commands.h"
#include "typedefs.h"

/*======function prototypes:=================================================*/

char		*ExtractToken_ (char *, int, char *, char *);
int		WriteLine_ (FILE *, RawAtomS *);
int		ExtractIndex_ (char *);
size_t		MainRefresh_ (MolComplexS *, int,
			      RuntimeS *, ConfigS *, GUIS *,
			      NearestAtomS *, size_t, unsigned int);
int		ControlRefresh_ (MolComplexS *, ConfigS *, GUIS *);

/*======execute group command:===============================================*/

int Group_ (MolComplexS *mol_complexSP, int mol_complexesN,
	    RuntimeS *runtimeSP, ConfigS *configSP, GUIS *guiSP,
	    NearestAtomS *nearest_atomSP, size_t pixelsN,
	    unsigned int *refreshIP, char *stringP)
{
size_t		raw_atom_struct_size;
int		max;
char		*remainderP;
char		tokenA[SHORTSTRINGSIZE];
FILE		*fileP;
int		chainI = 0;
static char	chainIDA[30] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
char		chainID;
int		complexID;
int		mol_complexI;
MolComplexS	*curr_mol_complexSP;
int		curr_mol_complexID;
size_t		atomI;
AtomS		*curr_atomSP;
RawAtomS	*raw_atomSP;
RawAtomS	copy_raw_atomS;
int		new_atomI = 0;
int		residue_sequenceI, old_residue_sequenceI = 0;
char		insertion_code, old_insertion_code = '\0';
int		new_residueI = 0;
int		arrayI = 0;
int		mol_complexIDA[MAX_GROUP_SIZE];
int		structures_in_groupN;
int		job_doneF = 0;
double		center_x = 0.0, center_y = 0.0, center_z = 0.0;
size_t		total_atomsN = 0;
double		reciprocal_total_atomsN;

/* The size of RawAtomS structure: */
raw_atom_struct_size = sizeof (RawAtomS);

/* Check the input string, there should be something: */
max = SHORTSTRINGSIZE;
if ((remainderP = ExtractToken_ (tokenA, max, stringP, " \t,;")) == NULL)
	{
	strcpy (runtimeSP->messageA,
		"Some additional parameters were expected!");
	runtimeSP->message_length = strlen (runtimeSP->messageA);
	return ERROR_NO_ID;
	}

/*---------------------------------------------------------------------------*/

/* Check for keyword SAVE: */
if (strstr (tokenA, "SAV") == tokenA)
	{
	/* If there is no group at all, there is nothing to save: */
	if (runtimeSP->groupF == 0)
		{
		strcpy (runtimeSP->messageA, "The group is not defined!");
		runtimeSP->message_length = strlen (runtimeSP->messageA);
		return ERROR_GROUP;
		}

	/* Now parse the original command line because the */
	/* copy of this line  was converted  to uppercase. */

	/* Skip the first and the second token: */
	remainderP = ExtractToken_ (tokenA, STRINGSIZE,
				    runtimeSP->curr_commandA, " \t\n");
	if (!remainderP) return ERROR_GROUP;
	remainderP = ExtractToken_ (tokenA, STRINGSIZE, remainderP, " \t\n");
	if (!remainderP) return ERROR_GROUP;

	/* The second token should be the name of the output file: */
	remainderP = ExtractToken_ (tokenA, STRINGSIZE, remainderP, " \t\n");
	if (!remainderP)
		{
		strcpy (runtimeSP->messageA, "Missing the output file name!");
		runtimeSP->message_length = strlen (runtimeSP->messageA);
		return ERROR_NO_FILE_NAME;
		}

	/* Try to open the output file: */
	if ((fileP = fopen (tokenA, "w")) == NULL)
		{
		strcpy (runtimeSP->messageA, "Failed to open output file!");
		runtimeSP->message_length = strlen (runtimeSP->messageA);
		return ERROR_OPEN_FAILURE;
		}

	/* If this point is reached, save the entire group to the same */
	/* PDB file.  The output file will not contain any header. The */
	/* structures forming  the group  will be  treated  as chains. */

	/* Scan each macromolecular complex: */
	for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
		{
		/* Prepare the pointer to the current macromol. complex: */
		curr_mol_complexSP = mol_complexSP + mol_complexI;

		/* Check the number of atoms;  for bad */
		/* and discarded complexes it is zero: */
		if (curr_mol_complexSP->atomsN == 0) continue;

		/* Prepare the chain identifier: */
		chainI %= 26;
		chainID = chainIDA[chainI];

		/* Increment the chain index: */
		chainI++;

		/* Atomic loop: */
		for (atomI = 0; atomI < curr_mol_complexSP->atomsN; atomI++)
			{
			/* Pointer to the current atom: */
			curr_atomSP = curr_mol_complexSP->atomSP + atomI;

			/* Pointer to raw atomic data: */
			raw_atomSP = &curr_atomSP->raw_atomS;

			/* Copy the raw atomic data: */
			memcpy (&copy_raw_atomS, raw_atomSP,
				raw_atom_struct_size);

			/* Replace the chain identifier: */
			copy_raw_atomS.chainID = chainID;

			/* Replace the atomic serial index: */
			copy_raw_atomS.serialI = new_atomI + 1;

			/* Update the new atomic index: */
			new_atomI++;

			/* Now take care about the residue serial number. */

			/* Copy the residue data associated with this atom: */
			residue_sequenceI = raw_atomSP->residue_sequenceI;
			insertion_code = raw_atomSP->residue_insertion_code;

			/* If this is the first atom  from this */
			/* structure, copy the residue sequence */
			/* index  and  residue  insertion code: */
			if (atomI == 0)
				{
				old_residue_sequenceI = residue_sequenceI;
				old_insertion_code = insertion_code;
				}

			/* If this atom  does not belong to the */
			/* same  residue as  the previous atom, */
			/* update the new residue serial index: */
			if ((residue_sequenceI != old_residue_sequenceI) ||
			    (insertion_code != old_insertion_code))
				{
				/* Update index: */
				new_residueI++;

				/* Update the residue sequence */
				/* index  and  insertion code: */
				old_residue_sequenceI = residue_sequenceI;
				old_insertion_code = insertion_code;
				}

			/* Replace the residue index: */
			copy_raw_atomS.residue_sequenceI = new_residueI + 1;

			/* Replace the residue insertion code: */
			copy_raw_atomS.residue_insertion_code = ' ';

			/* Write the output line: */
			WriteLine_ (fileP, &copy_raw_atomS);
			}
		}

	/* Close the output file: */
	fclose (fileP);

	/* Return to the caller, redrawing is not required: */
	return COMMAND_GROUP;
	}

/*---------------------------------------------------------------------------*/

/* Check for keyword SSC: */
else if (strstr (tokenA, "SSC") == tokenA)
	{
	/* If there is no group at all, there is nothing to save: */
	if (runtimeSP->groupF == 0)
		{
		strcpy (runtimeSP->messageA, "The group is not defined!");
		runtimeSP->message_length = strlen (runtimeSP->messageA);
		return ERROR_GROUP;
		}

	/* Now parse the original command line because the */
	/* copy of this line  was converted  to uppercase. */

	/* Skip the first and the second token: */
	remainderP = ExtractToken_ (tokenA, STRINGSIZE,
				    runtimeSP->curr_commandA, " \t\n");
	if (!remainderP) return ERROR_GROUP;
	remainderP = ExtractToken_ (tokenA, STRINGSIZE, remainderP, " \t\n");
	if (!remainderP) return ERROR_GROUP;

	/* The second token should be the name of the output file: */
	remainderP = ExtractToken_ (tokenA, STRINGSIZE, remainderP, " \t\n");
	if (!remainderP)
		{
		strcpy (runtimeSP->messageA, "Missing the output file name!");
		runtimeSP->message_length = strlen (runtimeSP->messageA);
		return ERROR_NO_FILE_NAME;
		}

	/* Try to open the output file: */
	if ((fileP = fopen (tokenA, "w")) == NULL)
		{
		strcpy (runtimeSP->messageA, "Failed to open output file!");
		runtimeSP->message_length = strlen (runtimeSP->messageA);
		return ERROR_OPEN_FAILURE;
		}

	/* If this point is reached, save the entire group to the same */
	/* PDB file.  The output file will not contain any header. The */
	/* structures forming  the group  will be  treated  as chains. */

	/* Scan each macromolecular complex: */
	for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
		{
		/* Prepare the pointer to the current macromol. complex: */
		curr_mol_complexSP = mol_complexSP + mol_complexI;

		/* Check the number of atoms;  for bad */
		/* and discarded complexes it is zero: */
		if (curr_mol_complexSP->atomsN == 0) continue;

		/* Prepare the chain identifier (space, for all structures): */
		chainID = ' ';

		/* Atomic loop: */
		for (atomI = 0; atomI < curr_mol_complexSP->atomsN; atomI++)
			{
			/* Pointer to the current atom: */
			curr_atomSP = curr_mol_complexSP->atomSP + atomI;

			/* Pointer to raw atomic data: */
			raw_atomSP = &curr_atomSP->raw_atomS;

			/* Copy the raw atomic data: */
			memcpy (&copy_raw_atomS, raw_atomSP,
				raw_atom_struct_size);

			/* Replace the chain identifier: */
			copy_raw_atomS.chainID = chainID;

			/* Replace the atomic serial index: */
			copy_raw_atomS.serialI = new_atomI + 1;

			/* Update the new atomic index: */
			new_atomI++;

			/* Now take care about the residue serial number. */

			/* Copy the residue data associated with this atom: */
			residue_sequenceI = raw_atomSP->residue_sequenceI;
			insertion_code = raw_atomSP->residue_insertion_code;

			/* If this is the first atom  from this */
			/* structure, copy the residue sequence */
			/* index  and  residue  insertion code: */
			if (atomI == 0)
				{
				old_residue_sequenceI = residue_sequenceI;
				old_insertion_code = insertion_code;
				}

			/* If this atom  does not belong to the */
			/* same  residue as  the previous atom, */
			/* update the new residue serial index: */
			if ((residue_sequenceI != old_residue_sequenceI) ||
			    (insertion_code != old_insertion_code))
				{
				new_residueI++;
				}

			/* Replace the residue index: */
			copy_raw_atomS.residue_sequenceI = new_residueI + 1;

			/* Replace the residue insertion code: */
			copy_raw_atomS.residue_insertion_code = ' ';

			/* Write the output line: */
			WriteLine_ (fileP, &copy_raw_atomS);
			}
		}

	/* Close the output file: */
	fclose (fileP);

	/* Return to the caller, redrawing is not required: */
	return COMMAND_GROUP;
	}

/*---------------------------------------------------------------------------*/

/* Check for keyword ALL: */
else if (strstr (tokenA, "ALL") == tokenA)
	{
	/* Set the catch flag to one for each macromolecular complex: */
	for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
		{
		/* Prepare the pointer to the current macromol. complex: */
		curr_mol_complexSP = mol_complexSP + mol_complexI;

		/* Check the number of atoms;  for bad */
		/* and discarded complexes it is zero: */
		if (curr_mol_complexSP->atomsN == 0) continue;

		/* Set the catch flag to one: */
		curr_mol_complexSP->catchF = 1;

		/* Mark this macromolecular complex as group member: */
		curr_mol_complexSP->group_memberF = 1;
		}

	/* Each macromolecular complex retains its own rotation center! */

	/* Return to the caller, redrawing is not required: */
	return COMMAND_GROUP;
	}

/*---------------------------------------------------------------------------*/

/* If this point is reached, keywords SAVE, SSC and ALL were */
/* not found.  Check are there  any identifiers  (integers). */

/* Parse the list of macromolecular identifiers: */
max = SHORTSTRINGSIZE;
remainderP = stringP;
while ((remainderP = ExtractToken_ (tokenA, max, remainderP, " \t,;")) != NULL)
	{
	/* Extract the macromolecular complex identifier */
	/* (an integer  larger than  zero  is expected): */
	if ((complexID = ExtractIndex_ (tokenA)) <= 0)
		{
		strcpy (runtimeSP->messageA,
			"Failed to extract the macromolecular complex");
		strcat (runtimeSP->messageA, " identifier!");
		runtimeSP->message_length = strlen (runtimeSP->messageA);
		return ERROR_NO_ID;
		}

	/* Store the macromolecular complex identifier: */
	mol_complexIDA[arrayI] = complexID;

	/* Increment and check the array index: */
	arrayI++;
	if (arrayI >= MAX_GROUP_SIZE) break;
	}

/* Copy the number of structures in a group: */
structures_in_groupN = arrayI;

/* If no identifiers were found, return: */
if (structures_in_groupN <= 0)
	{
	strcpy (runtimeSP->messageA,
		"Some additional parameters were expected!");
	runtimeSP->message_length = strlen (runtimeSP->messageA);
	return ERROR_NO_ID;
	}

/* If at least one identifier was found, reset all catch flags */
/* to zero  (this will  be changed later  for some complexes): */
for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
	{
	/* Prepare the pointer to the current macromolecular complex: */
	curr_mol_complexSP = mol_complexSP + mol_complexI;

	/* For all other complexes, set the catch flag to zero: */
	curr_mol_complexSP->catchF = 0;
	}

/* Set the catch flag for macromolecular */
/* complexes which belong  to the group: */
for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
	{
	/* Prepare the pointer to the current macromolecular complex: */
	curr_mol_complexSP = mol_complexSP + mol_complexI;

	/* Check the number of atoms;  for bad */
	/* and discarded complexes it is zero: */
	if (curr_mol_complexSP->atomsN == 0) continue;

	/* Copy the current macromolecular complex identifier: */
	curr_mol_complexID = curr_mol_complexSP->mol_complexID;

	/* Scan the group: */
	for (arrayI = 0; arrayI < structures_in_groupN; arrayI++)
		{
		/* If the current macromolecular complex is recognized: */
		if (curr_mol_complexID == mol_complexIDA[arrayI])
			{
			/* Set the catch flag to one: */
			curr_mol_complexSP->catchF = 1;

			/* Mark this macromolecular complex as group member: */
			curr_mol_complexSP->group_memberF = 1;

			/* Set the flag which says that group is */
			/* formed and break from the inner loop: */
			job_doneF = 1;
			break;
			}
		}
	}

/* If no complex was found, return negative value: */
if (job_doneF == 0)
	{
	strcpy (runtimeSP->messageA,
		"All complex identifiers are bad!");
	runtimeSP->message_length = strlen (runtimeSP->messageA);
	return ERROR_BAD_ID;
	}

/* Calculate the position of  the rotation center.  This rotation */
/* center will be common for all structures which form the group. */
/* The geometric center will be  calculated for all atoms and for */
/* all complexes,  including  the atoms  which are  not selected. */

/* Scan each macromolecular complex and take each atom into account: */
for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
	{
	/* Prepare the pointer to the current macromolecular complex: */
	curr_mol_complexSP = mol_complexSP + mol_complexI;

	/* Check the number of atoms;  for bad */
	/* and discarded complexes it is zero: */
	if (curr_mol_complexSP->atomsN == 0) continue;

	/* Atomic loop: */
	for (atomI = 0; atomI < curr_mol_complexSP->atomsN; atomI++)
		{
		/* Pointer to the current atom: */
		curr_atomSP = curr_mol_complexSP->atomSP + atomI;

		/* Add x, y and z to the totals: */
		center_x += curr_atomSP->raw_atomS.x[0];
		center_y += curr_atomSP->raw_atomS.y;
		center_z += curr_atomSP->raw_atomS.z[0];

		/* Increment the total number of atoms: */
		total_atomsN++;
		}
	}

/* Check the total number of atoms: */
if (total_atomsN == 0)
	{
	strcpy (runtimeSP->messageA, "There are no atoms in this group!");
	runtimeSP->message_length = strlen (runtimeSP->messageA);
	return ERROR_GROUP;
	}

/* Calculate the coordinates of the rotation center: */
reciprocal_total_atomsN = 1.0 / (double) total_atomsN;
center_x *= reciprocal_total_atomsN;
center_y *= reciprocal_total_atomsN;
center_z *= reciprocal_total_atomsN;

/* Replace the original rotation center vector with the group's center: */
for (mol_complexI = 0; mol_complexI < mol_complexesN; mol_complexI++)
	{
	/* Prepare the pointer to the current macromolecular complex: */
	curr_mol_complexSP = mol_complexSP + mol_complexI;

	/* Check is this complex a group member: */
	if (curr_mol_complexSP->group_memberF == 0) continue;

	/* If this point is reached, the current complex is a group member. */

	/* Backup the original rotation center vector: */
	curr_mol_complexSP->backup_vectorS.x =
			curr_mol_complexSP->rotation_center_vectorS.x;
	curr_mol_complexSP->backup_vectorS.y =
			curr_mol_complexSP->rotation_center_vectorS.y;
	curr_mol_complexSP->backup_vectorS.z =
			curr_mol_complexSP->rotation_center_vectorS.z;

	/* Replace the original coordinates: */
	curr_mol_complexSP->rotation_center_vectorS.x = center_x;
	curr_mol_complexSP->rotation_center_vectorS.y = center_y;
	curr_mol_complexSP->rotation_center_vectorS.z = center_z;
	}

/* Set the group flag in the RuntimeS structure: */
runtimeSP->groupF = 1;

/* Default editing mode (0, no editing) does not */
/* require  redrawing.  Other modes  require it: */
if (runtimeSP->edit_modeI != 0)
	{
	/* Reset the edit mode index: */
	runtimeSP->edit_modeI = 0;

	/* Refresh the main window: */
	(*refreshIP)++;
	MainRefresh_ (mol_complexSP, mol_complexesN,
		      runtimeSP, configSP, guiSP,
		      nearest_atomSP, pixelsN, *refreshIP);

	/* Refresh the control window: */
	ControlRefresh_ (mol_complexSP + runtimeSP->default_complexI,
			 configSP, guiSP);
	}

/* Return positive value on success: */
return COMMAND_GROUP;
}

/*===========================================================================*/