#!/usr/bin/env perl
use strict;
use lib $ENV{ 'MAIN_CONF'}  ;
use main_config ;
my %commands = () ;

if ( @ARGV <= 0 ) { push ( @{$commands{'help'}}, ' ' ) ; } 
else { (%commands) = main_config::get_commands ( @ARGV ) ; }

######################################################

my $MACRO_FILE = "center_next.com" ;
my $VALIDATE_FILE = "validate.dat" ;
my $TOP_LIST = "-1" ;
my $TASK = "NONE" ;
my $CENTER_BY = "OTHER" ;
my $CUT_OFF = 0 ;
my $FIELD = "AVERAGE" ;
my $FORM = "ANALYSIS" ;
my $PROC_ERROR = "FALSE" ;
my $doit = "FALSE" ;

my $keyword ; 
foreach $keyword ( sort keys (%commands) ) 
{
    if ( $keyword eq 'h' || substr($keyword, 0,4) eq "help" ) { goto HELP }
    if ( $keyword eq 'o' || substr($keyword, 0,6) eq "output" ) { goto OUTPUT }
    if ( $keyword eq 'i' || substr($keyword, 0,5) eq "input" ) { goto INPUT }
    if ( $keyword eq 'b' || substr($keyword, 0,2) eq "by" ) { goto CENTER_BY }
    if ( $keyword eq 'c' || substr($keyword, 0,6) eq "cutoff" ) { goto CUT_OFF }
    if ( $keyword eq 'f' || substr($keyword, 0,5) eq "field" ) { goto FIELD }
    if ( $keyword eq 'r' || substr($keyword, 0,4) eq "form" ) { goto FORM }
    if ( $keyword eq 'k' || substr($keyword, 0,4) eq "task" ) { goto TASK }
    if ( $keyword eq 't' || substr($keyword, 0,3) eq "top" ) { goto TOP }
    if ( substr($keyword, 0,4) eq "doit" ) { $doit = "TRUE" ; next ; }
    print "unknown option >$keyword<\n" ; $PROC_ERROR = "TRUE" ;

HELP:
    print << "END_OF_TEXT" ;
  -h|--help)    prints this message with available options and current status 
  -o|--output)  set file name of the created macro [$MACRO_FILE]
  -i|--input)   set file name of the validation input [$VALIDATE_FILE]
  -c|--cutoff)  cutoff value [$CUT_OFF]
  -f|--field)   field [AVERAGE/RMS/MIN/MAX] [$FIELD]
  -r|--form)    field [ANALYSIS/PEAKS/PDB/POINT] [$FORM]
  -k|--task)    task [NONE/MINIMIZE/FIT_SIDE/AUTO_FIT/PICK_ATOM/MY_TASK] [$TASK]
  -t|--top)     top [$TOP_LIST]
  -b|--by)      task [GEO_BOND/GEO_ANGLE/PACK/ENE_DENS/OTHER] [$CENTER_BY]
     --doit)    create the macro
END_OF_TEXT
    exit;

OUTPUT: $MACRO_FILE = shift ( @{$commands{$keyword}} ) ; next ;
INPUT: $VALIDATE_FILE = shift ( @{$commands{$keyword}} ) ; next ;
CUT_OFF: $CUT_OFF = shift ( @{$commands{$keyword}} ) ; next ;
FIELD: $FIELD = shift ( @{$commands{$keyword}} ) ; next ;
FORM: $FORM = shift ( @{$commands{$keyword}} ) ; next ;
TOP: $TOP_LIST = shift ( @{$commands{$keyword}} ) ; next ;
TASK: $TASK = shift ( @{$commands{$keyword}} ) ; next ;
CENTER_BY: $CENTER_BY = shift ( @{$commands{$keyword}} ) ; next ;
}
if ( $PROC_ERROR eq "TRUE" ) { exit ; }

#######################################

if ( $doit eq "FALSE" )  { exit ; }

######################################

if ( $FORM eq "ANALYSIS" ) { create_analysis_list () ; }
if ( $FORM eq "PEAKS" ) { create_peak_list () ; }
if ( $FORM eq "PDB" ) { create_pdb_list () ; }

exit ;

sub create_pdb_list
{

    my $PDB_FILE = $VALIDATE_FILE ;

    my ( $ref_segm_list, $ref_resid_list, $ref_seq_list, $ref_unit_cell, $ref_space_group, $ref_atom_list, $ref_xyz_list, $ref_b_list, $ref_w_list, $ref_atom_el_list, $ref_helix_list, $ref_beta_list ) 
	= main_config::read_pdb_file ( $PDB_FILE ) ;
#my ( $ref_segm_list, $ref_resid_list, $ref_seq_list, $ref_unit_cell, $ref_space_group, $ref_atom_list, $ref_xyz_list, $ref_b_list, $ref_w_list, $re#f_atom_el_list ) = 
#    main_config::read_pdb_file ( $PDB_FILE ) ;

    my @my_segm_list = @{$ref_segm_list} ;
    my @my_resid_list = @{$ref_resid_list} ;
    my @my_seq = @{$ref_seq_list} ;
    my @my_xyz_list = @{$ref_xyz_list} ;
    my @my_w_list = @{$ref_w_list} ;
    my @my_b_list = @{$ref_b_list} ;

    open ( MACRO,  "> $MACRO_FILE" ) || die "FILE OPEN ERROR: >$MACRO_FILE< \n" ;

#    my %atom_list = () ;

#    my $xx ;
#    for ( $xx = 0; $xx < scalar(  @my_b_list  ) ; $xx++ )
#    {
#	 push( @{$atom_list{$my_b_list[$xx]}}, $xx );
#    }
#
#    my @key_list = () ;
#    my $key ;
#    foreach $key ( keys (%atom_list) ) 
#    { 
#	push(@key_list, $key ) ; 
##	print $key, " = > ", join(" ", @{$atom_list{$key}}), "\n" ;
#    }


#    my @sorted = sort { $b <=> $a } @key_list ;
    my $count = 0 ;


##    foreach $key (@sorted)  
#    {
##foreach $key ( sort numerically keys (%aver_list) ) 
##    print " key:", $key, " ++ ", join(" ", @{$atom_list{$key}}), "\n" ;
#	my $LIST ;
#	foreach $LIST ( @{$atom_list{$key}} )
#	{
#	    $count++ ;
#    printf MACRO "! %d %d weight %9.3f \n", 
#	    $count, $key, $my_b_list[$LIST] ;
#    printf MACRO "image center coor %9.3f %9.3f %9.3f \n", 
#	    $my_xyz_list[$LIST*3], 
#	    $my_xyz_list[$LIST*3+1], $my_xyz_list[$LIST*3+2] ;
#    print MACRO "image map over \n" ;
#    print MACRO "return \n\n" ;
#	}
##	if ( $TOP_LIST > 0 && $count >= $TOP_LIST ) { last ; }
#	if ( $TOP_LIST > 0 && $count >= $TOP_LIST ) { 

    my $LIST=0 ;
    while ( $LIST < $#my_xyz_list )
    {
	$count++ ;
	printf MACRO "! %d \n", $count ;
printf MACRO "image center coor %9.3f %9.3f %9.3f \n", 
	$my_xyz_list[$LIST], $my_xyz_list[$LIST+1], $my_xyz_list[$LIST+2] ;
print MACRO "image hist center \n" ;
print MACRO "image map over \n" ;
print MACRO "return \n\n" ;
	$LIST = $LIST+3 ;
    }
#    if ( $TOP_LIST > 0 && $count >= $TOP_LIST ) { last ; }
#    if ( $TOP_LIST > 0 && $count >= $TOP_LIST ) 
#   { 
	if ( $count > 0 ) 
	{
print MACRO << "END_OF_TEXT" ; 
echo \" CENTER_NEXT> END OF LIST REACHED - START FROM THE TOP\"
rewind unit 66
END_OF_TEXT
	}
close ( MACRO ) ;
}

sub create_peak_list
{

    my $PDB_FILE = "atom_center_list.pdb" ;

    my ( $ref_segm_list, $ref_resid_list, $ref_seq_list, $ref_unit_cell, $ref_space_group, $ref_atom_list, $ref_xyz_list, $ref_b_list, $ref_w_list, $ref_atom_el_list, $ref_helix_list, $ref_beta_list ) 
	= main_config::read_pdb_file ( $PDB_FILE ) ;
#    my ( $ref_segm_list, $ref_resid_list, $ref_seq_list, $ref_unit_cell, $ref_space_group, $ref_xyz_list, $ref_w_list, $ref_b_list ) = 
#	main_config::read_pdb_file ( $PDB_FILE ) ;

    my @my_segm_list = @{$ref_segm_list} ;
    my @my_resid_list = @{$ref_resid_list} ;
    my @my_seq = @{$ref_seq_list} ;
    my @my_xyz_list = @{$ref_xyz_list} ;
    my @my_w_list = @{$ref_w_list} ;
    my @my_b_list = @{$ref_b_list} ;

    my %atom_list = () ;

    my $xx ;
    for ( $xx = 0; $xx < scalar(  @my_b_list  ) ; $xx++ )
    {
	 push( @{$atom_list{$my_b_list[$xx]}}, $xx );
    }

    my @key_list = () ;
    my $key ;
    foreach $key ( keys (%atom_list) ) 
    { 
	push(@key_list, $key ) ; 
    }

open ( MACRO,  "> $MACRO_FILE" ) || die "FILE OPEN ERROR: >$MACRO_FILE< \n" ;

    my @sorted = sort { $b <=> $a } @key_list ;
    my $count = 0 ;
    foreach $key (@sorted)  
    {
#foreach $key ( sort numerically keys (%aver_list) ) 
#    print " key:", $key, " ++ ", $aver_list{$key}, "\n" ;
	my $LIST ;
	foreach $LIST ( @{$atom_list{$key}} )
	{
	    $count++ ;
    printf MACRO "! %d %d weight %9.3f \n", 
	    $count, $key, $my_b_list[$LIST] ;
    printf MACRO "image center coor %9.3f %9.3f %9.3f \n", 
	    $my_xyz_list[$LIST*3], 
	    $my_xyz_list[$LIST*3+1], $my_xyz_list[$LIST*3+2] ;
    print MACRO "image hist center \n" ;
    print MACRO "image map over \n" ;
    print MACRO "return \n\n" ;
	}
	if ( $TOP_LIST > 0 && $count >= $TOP_LIST ) { last ; }
    }
    
if ( $count > 0 ) 
{
print MACRO << "END_OF_TEXT" ; 
echo \" CENTER_NEXT> END OF LIST REACHED - START FROM THE TOP\"
rewind unit 66
END_OF_TEXT
}
close ( MACRO ) ;
}

sub create_analysis_list
{

open ( TXT,  "< $VALIDATE_FILE" ) || die "FILE OPEN ERROR: >$VALIDATE_FILE< \n" ;
my %field_name = () ;
if ( $CENTER_BY eq "GEO_BOND" ) 
{
    $field_name{"RMS"} = 2 ;
    $field_name{"MIN"} = 4 ;
    $field_name{"MAX"} = 6 ;
    $field_name{"AVERAGE"} = 8 ;
}
else
{
    $field_name{"AVERAGE"} = 1 ;
    $field_name{"RMS"} = 3 ;
    $field_name{"MIN"} = 5 ;
    $field_name{"MAX"} = 7 ;
}

my $field = $field_name{$FIELD} ;
my %aver_list = () ;
my $key ;
my @words = () ;

my @file_lines = () ;

while(<TXT>)
{
    chomp ;
    my $line = $_ ;
    push( @file_lines, $line ) ;
    @words = split(" ", $line) ;
    $key = $words[$field] ;
    if ( $CENTER_BY eq "GEO_BOND" ) { $key = abs($key) ; }
    
#    print " my key $key =+= $field \n" ;
#    print $words[0], " read in :", $field, ", ", $words[9], ", ", $words[8], ", ", $words[7],", ", $words[1], "\n" ;

    if ( $words[0] > 0 ) 
    {
	if ( $key =~ m"\*" ) { $key = 99999999 ; }
	if ( substr($key, 0, 3) eq "..." ) { $key = -99999999 ; }
	if ( $CENTER_BY eq "GEO_BOND" )
	{
	    push( @{$aver_list{$key}}, $words[10] );
	}
	else
	{
	    push( @{$aver_list{$key}}, $words[9] );
	}
    }
}

print "creating \"$MACRO_FILE\" \n" ;
if ( -f $MACRO_FILE ) 
{ 
    my $back = join ('', $MACRO_FILE, ".old" ) ; 
    system ( "mv  $MACRO_FILE $back" ) ; 
}

open ( MACRO,  "> $MACRO_FILE" ) || die "FILE OPEN ERROR: >$MACRO_FILE< \n" ;

my ( $summ, $summ2, $nn ) ;
my $rmsa=0 ; my $aver=0 ;
my $tollerance ;
my @key_list = () ;
$summ = 0; $summ2 = 0 ; $nn = 0 ; $tollerance = -9999 ;
foreach $key ( keys (%aver_list) ) 
{ 
    push(@key_list, $key ) ; 
    $summ = $summ+$key; 
    $summ2 = $summ2 + $key*$key ;
    $nn++ ;
}
$rmsa = 0 ; $aver = 0 ;
if ( $nn > 1 ) 
{
    $rmsa = abs($summ2 - $summ / $nn * $summ ) ;
    $rmsa = sqrt($rmsa / ( $nn-1) ) ;
    $aver = $summ / $nn ;
    if ( $CUT_OFF != 0 ) 
    { 
	$tollerance = $summ / $nn + $rmsa * $CUT_OFF ; 
#print " tollerance:", $tollerance, "\n" ;
    }
    if ( $CENTER_BY eq "GEO_BOND" || $CENTER_BY eq "ANGLE" || 
	 $CENTER_BY =~ "PACK"  || $CENTER_BY eq "ENE_DENS" ) 
    { 
	$rmsa = 1 ; 
	$aver = 0 ; 
	$tollerance = $CUT_OFF ;
    }
}
else
{
    $rmsa = 1 ; 
    print " NN == $nn \n" ;
}
if ( $CENTER_BY ne "NONE" ) { $tollerance = $rmsa * $CUT_OFF; }


if ( $FIELD eq "MIN" && $CENTER_BY eq "ENE_DENS" )
{ 
#    $tollerance = - abs($tollerance) ; 
}
print " tolerance: $tollerance \n" ;


my @sorted = () ;
if ( $FIELD eq "MIN" )
{ 
    print " reverse sort:\n" ;
    @sorted = reverse sort { $b <=> $a } @key_list ; 
}
else 
{ 
    print " normal sort:\n" ;
    @sorted = sort { $b <=> $a } @key_list ; 
}

my $count = 0 ; my $done = 0 ;
foreach $key (@sorted)  
{
#foreach $key ( sort numerically keys (%aver_list) ) 
#    print " key:", $key, " ++ ", join(" ", @{$aver_list{$key}}), "\n" ;
    my $LIST ;
    foreach $LIST ( @{$aver_list{$key}} )
    {
	$count++ ;
#	print " tolerance:", $tollerance, " - ", $key, "\n" ;
#	if ( ( $TOP_LIST < 0 || $count <= $TOP_LIST ) &&
#	     ( $tollerance < -9990 || $key >= $tollerance ) )
	my $ACCEPT_THIS_ONE = "NOP" ;
	if ( $TOP_LIST < 0 || $count <= $TOP_LIST ) 
	{
	    if ( $FIELD eq "MIN" )
	    { 
		if ( $tollerance > -9990 && $key <= $tollerance ) 
		{ # print "MIN: $key\n"; 
		  $ACCEPT_THIS_ONE = "YES" ; }
	    }
	    else
	    {
		if ( $tollerance > -9990 && $key >= $tollerance ) 
		{ # print "MAX: $key\n"; 
		  $ACCEPT_THIS_ONE = "YES" ; }
	    }
	}

if ( $ACCEPT_THIS_ONE eq "YES" )
{
#    print " count $count top $TOP_LIST $key \n" ;
    $done++ ;
    my  $tmp = ( $key - $aver ) ; 
    my  $tmp_sig = $tmp / $rmsa ;
printf MACRO "\n! %9.3f aver %9.3f delta %9.3f sigma %9.3f delta/sigma %9.3f",  
    $key, $aver, $tmp, $rmsa, $tmp_sig ;

print MACRO << "END_OF_TEXT" ; 

set vari NEXT_RESI int global = $LIST
set vari II int atom sele history .and. atom numb NEXT_RESI end
if ( II .ne. NEXT_RESI ) image hist pick NEXT_RESI
END_OF_TEXT

if ( $CENTER_BY eq "ENE_VDW" || $CENTER_BY =~ "PACK" )
{ print MACRO "image hist neigh 3.0\n" ; }

if ( $TASK eq "MINIMIZE" ) 
{
print MACRO << "END_OF_TEXT" ; 
key active sele by resi atom numb NEXT_RESI end
key passive sele by bond  active  end
minimize sele active end sele passive end
minimize sele active end sele passive end
minimize sele active end sele passive end
END_OF_TEXT
}

if ( $TASK eq "FIT_SIDE" ) 
{
print MACRO << "END_OF_TEXT" ; 

key active sele by resi atom numb NEXT_RESI end
key passive sele by bond  active  end

< MAIN_CMDS_NEW:auto_stuff_side.cmds WORK_SEGM FIT
minimize sele active end sele passive end
minimize sele active end sele passive end
END_OF_TEXT
}

if ( $TASK eq "AUTO_FIT" ) 
{
print MACRO << "END_OF_TEXT" ; 

key active sele by resi atom numb NEXT_RESI end
key passive sele by bond  active  end

< MAIN_CMDS_NEW:auto_fit.cmds 
minimize sele active end sele passive end
minimize sele active end sele passive end
END_OF_TEXT
}

if ( $TASK eq "FIT_PEPT" ) 
{
print MACRO << "END_OF_TEXT" ; 

key active sele by resi atom numb NEXT_RESI end
key passive sele by bond  active  end

< MAIN_CMDS_NEW:auto_stuff_find_pept.cmds 
minimize sele active end sele passive end
minimize sele active end sele passive end
END_OF_TEXT
}

if ( $TASK eq "PICK_ATOM" ) 
{
print MACRO << "END_OF_TEXT" ; 

key active sele by resi atom numb NEXT_RESI end
key passive sele by bond  active  end
image hist pick NEXT_RESI

END_OF_TEXT
}

if ( $TASK eq "MY_TASK" ) 
{
print MACRO << "END_OF_TEXT" ; 

key active sele by resi atom numb NEXT_RESI end
key passive sele by bond  active  end

< my_task.com

END_OF_TEXT
}

print MACRO << "END_OF_TEXT" ; 
return
END_OF_TEXT

	}
    }
}
#print " MY NN: ++++++++++++++++= $nn $count $done\n" ;
if ( $done > 0 ) 
{
print MACRO << "END_OF_TEXT" ; 
echo \" CENTER_NEXT> END OF LIST REACHED - START FROM THE TOP\"
rewind unit 66
END_OF_TEXT
}


close ( MACRO ) ;
}
