#!/usr/bin/perl -w

use strict;

my $numArgs = $#ARGV + 1;  # calculate number of arguments passed to program from command line
my $programchoice = $ARGV[0];
our $inputfile;
our $outputfile;  # initialise filename variables
our @fpmagmax;
our @fpmagmin;
our $fpmagmax;
our $fpmagmin;
our @fpmagsum;
our $fpmag;
our $molindex = $ARGV[1];

for my $i (0 .. 3) {
	$fpmagmin[$i] = 0;
	$fpmagmax[$i] = 0; # easier way to do this?
}

if ($numArgs < 2) {
	die "Please use a program option and a molecule index in the arguments";
}

elsif ($numArgs >= 4) {
	$inputfile = "$ARGV[2]";
	print "Input filename received from command line.  Opening '$inputfile'...\n";
	
	$outputfile = "$ARGV[3]";
	print "Output filename received from command line.  Opening '$outputfile'...\n";
}

elsif ($numArgs == 3) {
	$inputfile = "$ARGV[2]";
	print "Input filename received from command line.  Opening '$inputfile'...\n";
	
	print "Please enter the path to your output file: ";
	chomp($outputfile = <STDIN>);
}

else {
	print "Please enter the path to your input file: ";
	chomp($inputfile = <STDIN>);
	
	print "Please enter the path to your output file: ";
	chomp($outputfile = <STDIN>);

}	

open INPUTFILE, "<$inputfile" || die "An error occured whilst opening your input file '$inputfile': $!";

if ($programchoice == 2 || $programchoice == 3) {
	open OUTPUTFILE, ">>$outputfile" || die "An error occured whilst creating your output file '$outputfile': $!";
}

else {
	open OUTPUTFILE, ">$outputfile" || die "An error occured whilst creating your output file '$outputfile': $!";
}

print "INPUT: $inputfile\n";
print "OUTPUT: $outputfile\n";  #display input and output file names

eval {
	if ($programchoice == 1) {

		&xedparse;
	}

	elsif ($programchoice == 2) {

		&conformerselect_min;
	}
	
	elsif ($programchoice == 3) {

		&conformerselect_boltzmann;
	}
	
	else {
		die "your program option is invalid";
	}
	
	
		
};

close INPUTFILE;
close OUTPUTFILE;  # close input and output files


sub xedparse {

	print "XEDparse\n";
	
	my @conformers;
	my @energy;
	my $energy;

	my @fields;
	my $fields;
	my $lineindex = 0;

	my $endofconf = 1;
	my $confindex = 0;
	
	while(defined(my $line = <INPUTFILE>)) {  # read through each line in turn while they still exist 
		
		$_ = $line;
		chomp $_;
		
		@fields = split;  # split line into fields according to whitespace delimiter
		
		if ($fields[5]) {
			$fpmag = sqrt($fields[5] * $fields[5]);
		}
		
		# repetitive code 
				
		if ($fields[0] == -5) {
			
			# line contains a negative field point
			
			&fpparse(0);
				
		} 
		
		elsif ($fields[0] == -6) {
			
			# line contains a positive field point
			
			&fpparse(1);
				
		} 
		
		elsif ($fields[0] == -7) {
			
			# line contains a van der Waals field point
			
			&fpparse(2);
				
		} 
		
		elsif ($fields[0] == -8) {
			
			# line contains a steric field point
			
			&fpparse(3);
				
		}
		
		elsif ($fields[0] == 0 && $fields[1] == 0 && $fields[2] == 0 && $fields[3] == 0) {
			
			# end of conformer so print the details to the output, begin a new line and reset variables
			print OUTPUTFILE "$confindex,";
			$endofconf = 1;
			$confindex++;
			
			for my $i (0 .. 3) {
				print OUTPUTFILE "$fpmagsum[$i],";
				$fpmagsum[$i] = 0;
				$fpmagmax[$i] = 0;
				$fpmagmin[$i] = 0;
			}
			
			print OUTPUTFILE "$energy\n";
		}
		
		elsif ($endofconf == 1) {
			# first line of new conformer so extract energy
			
			$endofconf = 0;
			$energy = $fields[0];
			
		}	
		
		$lineindex++;  # increment current line

	}
	
	return;
}

sub fpparse {
	
	if ($fpmagmin[$_[0]] == 0 || $fpmagmax[$_[0]] == 0) {
		$fpmagmin = $fpmag;
		$fpmagmax = $fpmag;
	}
			
	if ($fpmag < $fpmagmax[$_[0]]) {
		$fpmagmax[$_[0]] = $fpmag;
	}
			
	if ($fpmag > $fpmagmin[$_[0]]) {
		$fpmagmin[$_[0]] = $fpmag;
	}
			
	$fpmagsum[$_[0]] += $fpmag;
				
	return;
}


sub conformerselect_min {

	print "ConformerSelect - Minimum Energy\n";

	my @fields;
	my $fields;
	my $min_energy = 0;
	my $min_energy_lineindex = 0;
	my $lineindex = 0;
	
	if (eof INPUTFILE) {
		print "blank file!\n";
		print OUTPUTFILE ",,,,$molindex,\n"
	}
	
	while(defined(my $line = <INPUTFILE>)) {  # read through each line in turn while they still exist 
		
		$line =~ s/[\n]+//g;  # remove new line characters
		@fields = split(/,\t*/, $line);  # split line into fields according to comma delimiter
	
		if ($min_energy == 0) {
			$min_energy = $fields[5];
		}
	
		if ($fields[5] < $min_energy) {
			$min_energy = $fields[5];
			$min_energy_lineindex = $lineindex;
		}	
		
		$lineindex++;
	}
	
	

	
	$lineindex = 0;
	close INPUTFILE;
	open INPUTFILE, "<$inputfile" || die "An error occured whilst opening your input file '$inputfile': $!"; # reset input file
	
	#my $lines = 0;
	#while (defined(my $line = <INPUTFILE>)){
	#	$lines++;
	#}
	
	close INPUTFILE;
	open INPUTFILE, "<$inputfile" || die "An error occured whilst opening your input file '$inputfile': $!"; # reset input file
	
	#print OUTPUTFILE "Negative,Positive,VdW,Steric,Molecule Index,Conformer Index\n";
	

			
	while(defined(my $line = <INPUTFILE>)) {  # read through each line in turn while they still exist 
		
		$line =~ s/[\n]+//g;  # remove new line characters
		@fields = split(/,\t*/, $line);  # split line into fields according to comma delimiter
		#my $numberfields = @fields;
		#print "\$lines: $lines\n";
		#print "$numberfields\n";
	
		if ($lineindex == $min_energy_lineindex) {
		
			for my $i (1 .. 4) {	
				if (!$fields[$i]) {
					$fields[$i] = 0;
				}
				#print "$fields[$_]";
			}
		
			print OUTPUTFILE "$fields[1],$fields[2],$fields[3],$fields[4],$molindex,$lineindex\n";
		}	
		
		$lineindex++;
	}
	
	return;
}

sub conformerselect_boltzmann {

	print "ConformerSelect - Boltzmann Average\n";

	my @fields;
	my $fields;
	my $partition_element = 0;
	my $partition_function = 0;
	my $boltzmann = 0.0019871;
	my $temp = 1000;
	my @descriptor_average;
	my $descriptor_average;
	my $boltzmann_distfn;
	my $energy;
	
	while(defined(my $line = <INPUTFILE>)) {  # read through each line in turn while they still exist 
		
		$line =~ s/[\n]+//g;  # remove new line characters
		@fields = split(/,\t*/, $line);  # split line into fields according to comma delimiter
		$energy = $fields[5];
	
		$partition_element = exp(-$energy / ($boltzmann * $temp));
		$partition_function += $partition_element;

	}
	
	close INPUTFILE;
	open INPUTFILE, "<$inputfile" || die "An error occured whilst opening your input file '$inputfile': $!"; # reset input file
	
	#print OUTPUTFILE "Negative,Positive,VdW,Steric,Molecule Index\n";
		
	while(defined(my $line = <INPUTFILE>)) {  # read through each line in turn while they still exist 
		
		$line =~ s/[\n]+//g;  # remove new line characters
		@fields = split(/,\t*/, $line);  # split line into fields according to comma delimiter
		$energy = $fields[5];
		#print "$energy / ($boltzmann * $temp)\n";
		$partition_element = exp(-$energy / ($boltzmann * $temp));
		#print "$partition_element / $partition_function\n";
		$boltzmann_distfn = ($partition_element / $partition_function);

		for my $i (0 .. 3) {
			$descriptor_average[$i] += ($boltzmann_distfn * $fields[$i +1]);
		}
		
	}
	
	print OUTPUTFILE "$descriptor_average[0],$descriptor_average[1],$descriptor_average[2],$descriptor_average[3],$molindex\n";
	
	return;
}