#!/usr/local/bin/perl -w
# given pdb file uses A carbons to find closes amino acids to CA
# finds closes set that will enclose CA in its convex hull

$firstTime=1;
$linecount=0; 
open (INPUT, $ARGV[0])|| die "couldnt open input";
$_=$ARGV[0];
(/data\/pdb(.*)\.ent/)||die "didnt match";
open (OUTPUT, ">results/".$1.".localhull")|| die "couldnt open output";

# find a CA to center on and store it's x,y,z
$hitca=0;
while (<INPUT>){
    if (/HETATM\s+.+CA\s+CA\s+.*(\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
	#print "hit CA";
	$hitca=1;
	$CAx=$2;
	$CAy=$3;
	$CAz=$4;
    }
}
if (!$hitca) {die "couldnt find ca\n";}
print "CA location is ".$CAx." ".$CAy." ".$CAz."\n";
close INPUT;
open (INPUT, $ARGV[0])|| die "couldnt open input";

$index=0;

while (<INPUT>){
# scan atom's store distance name and x,y,z

   #if(/ATOM\s+.+\s+CA\s+(\S+)\s(.+)\s(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
   #if(/ATOM\s+.+\s+CA\s+(\S+)\s(\S+\s*\d*)\s(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
   #if(/ATOM\s+.+\s+CA\s+(\S+)\s(.+)\s(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
   # if(/ATOM\s+.+\s+CA\s+(\S+)\s.+\s(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
#if(/ATOM\s+\S+\s+\S+\s+(\S+)\s+(.+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)\s+.+/){
if(/ATOM\s+\S+\s+\CA\s+(\S+)\s+(.+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)\s+.+/){
       #print "value is ".$2." ".$3." ".$4."\n";
       $moleculename[$index]=$2;
       $xval[$index]=$3;
       $yval[$index]=$4;
       $zval[$index]=$5;
       $distance[$index]=($3-$CAx)**2+($4-$CAy)**2+($5-$CAz)**2;
       $distance[$index]=$distance[$index]**.5;

       # this has index in array for each distance
       $atomsHash{$distance[$index]}=$index;
       $name[$index]=$1;
       $index++;
       
	
    }
  
}
close INPUT;
# sort list of distances numerically we can then go over distance array using is as key for hash to get
# the index of the corresponding name and atom records
@sortedlist=sort { $a <=> $b } @distance;

$counter=1;
$done=0;

while (!$done){
    #write set to file
    initializealtfile($counter);

    #check if it encloses
    if (enclosed()){
	print "enclosed\n";
	printresults();
	$done=1;
    }

    else { 
	#check for case where we had error with qhull
	if ($done){
	    print OUTPUT "QHULL ERROR\n";
	}
	$counter++;
    }

    if ($counter>$#xval) {
	# write fact that doesnt enclose 
	print "final set doesn't enclose\n";
	print OUTPUT "DOESNT ENCLOSE\n";
	$done=1;
    }
}



close INPUT;
close OUTPUT;
# remove temp file for this run
unlink("not_localpoints_ca.qhull");

#end main


#prints names of all amino acids who's Acarbons were included in latest shell calc
sub printresults{
    for ($it=0;$it<$counter;$it++){
	$valindex=$atomsHash{$sortedlist[$it]};
	print OUTPUT $name[$valindex]."\n";
	print $name[$valindex]."\n";
    }
}


# checks if the current set of atoms forms a convex hull that encloses the CA ion

sub enclosed{
# run qhull on file if CA's in hull false if not return true
# Alot of times these points are close to coplanar the -QJ jiggles the points so theyre
# no longer so the -Pp supresses the warning message when the hull is narrow
# because the points were close to coplanar  

if ($linecount<4){ return 0;}

$test=system("cat localpoints_ca.qhull|./qhull/qhull/qhull -Fx -QJ -Pp>temp.txt");
if ($test!=0) {$done=1};

open (OUT1, "temp.txt")|| die "couldnt open qhull results";
#print $counter."\n";
while (<OUT1>){
   # print $_."\n";
    
    #chomp;
    # if index=counter is the CA if it appears in shell its not enclosed
    if($_==$counter){
		#print "false";
		return 0;}
}
return 1; 
}




# writes counter elements from sorted list to file plus the ca position
#uses sorted distance list in hash to get index for x,y,z values
# is a somewhat dangerous assumption that all distances are different
#if theyre not some amino acids dissapeare from consideration 
sub initializealtfile
{    
    open (POINTFILE, ">localpoints_ca.qhull");
   
    if ($firstTime){$linecount=1;} #inc for CA if first time
    

    $newstuffindex=0;

    @newlyaddedstuff=();

    #find new entries values
# dont need this since starting from zero 
    #if ($firstTime){
	
#	for ($x=0;$x<($_[0]-1);$x++){
#	     $valindex=$atomsHash{$sortedlist[$x]};
#	     addnewvaluestolist($valindex);
#	}
#    }
  
    $valindex=$atomsHash{$sortedlist[($_[0]-1)]};
    addnewvaluestolist($valindex);

    # print header
    print POINTFILE "3\n".$linecount."\n";
   

    # print old stuff to file
    if (!$firstTime){
	 #this one  holds the data without header lines so it can be added to 
	open (NOTPOINTFILE, "not_localpoints_ca.qhull")||die "cant open notpointifile";
	while (<NOTPOINTFILE>){
	    print POINTFILE $_;
	    }
	 close NOTPOINTFILE;
    }
   $firstTime=0;

    open (NOTPOINTFILE,">>not_localpoints_ca.qhull");


    #while (<@newlyaddedstuff>){
#	print POINTFILE $_;
#	print NOTPOINTFILE $_;
#    }

    for ($b=0;$b<@newlyaddedstuff;$b++){
	print POINTFILE $newlyaddedstuff[$b];
	print NOTPOINTFILE $newlyaddedstuff[$b];
	#print "added\n";
    }

    # add ca on last line depends on this 
    print POINTFILE $CAx." ".$CAy." ".$CAz."\n";
	
    close NOTPOINTFILE;
    close POINTFILE;

}

# given index takes its name searches for entries with that name and adds them to newlyadded data
sub addnewvaluestolist
{
    $thename=$moleculename[$_[0]];
    open (INPUT, $ARGV[0])|| die "couldnt open input in addnevaluestolist"; 

    while (<INPUT>){
	#if(/ATOM\s+\S+\s+\S+\s+(\S+)\s+(.+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+\.\d+\s+\d+\.\d+.*/){
	#if(/ATOM\s+.+\s+CA\s+(\S+)\s(.+)\s(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
	# if(/ATOM\s+.+\s+CA\s+(\S+)\s(\S+\s*\d*)\s(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
	#if(/ATOM\s+.+\s+CA\s+(\S+)\s(.+)\s(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+\d+.+\d+.+\d+.+/){
	if(/ATOM\s+\S+\s+\S+\s+(\S+)\s+(.+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)\s+.+/){
	#print "match ".$2."\n";
	    if ($2 eq $thename){
		$newlyaddedstuff[$newstuffindex]=$3." ".$4." ".$5."\n";
		$newstuffindex++;
		$linecount++;
		#print "added in addnevalue".$2."\n";
	    }
	}
    }
    close INPUT;
}
