#!/usr/bin/perl

if(@ARGV<2){
  printf("usage: z-score.pl sulfated_list nonsulfated_list [seed]\n");
  exit(1);
}

if(@ARGV>2){
  srand($ARGV[2]);
}

open(IN,"<$ARGV[0]");
my $index=0;
while(<IN>){
  chomp();
  my @tmp = split(" ",$_);
  if(substr($tmp[0],0,1) ne "#"){
    $positive[$index][0]=$tmp[2];
    $positive[$index][1]=$tmp[3];
    $positive[$index][2]=$tmp[4];
    $positive[$index][3]=$tmp[6];
    $positive[$index][4]=$tmp[7]/100;
    $index++;
  }
}
my $np = $index;
close(IN);

open(IN,"<$ARGV[1]");
$index=0;
while(<IN>){
  chomp();
  my @tmp = split(" ",$_);
  if(substr($tmp[0],0,1) ne "#"){
    $negative[$index][0]=$tmp[2];
    $negative[$index][1]=$tmp[3];
    $negative[$index][2]=$tmp[4];
    $negative[$index][3]=$tmp[6];
    $negative[$index][4]=$tmp[7]/100;
    $index++;
  }
}
my $nn = $index;
close(IN);

my @coef;
my $max_c = 60;

$coef[1]=1;
$coef[0]=$coef[2]=$coef[3]=$coef[4]=0;
&compute_z();
printf("HELIX: %f\n", $zscore);

$coef[2]=1;
$coef[0]=$coef[1]=$coef[3]=$coef[4]=0;
&compute_z();
printf("STRAND: %f\n", $zscore);

$coef[3]=1;
$coef[0]=$coef[1]=$coef[2]=$coef[4]=0;
&compute_z();
printf("COIL: %f\n", $zscore);

$coef[4]=1;
$coef[0]=$coef[1]=$coef[2]=$coef[3]=0;
&compute_z();
printf("EXPOSURE: %f\n", $zscore);


$coef[0]=1;
$coef[1]=$coef[2]=$coef[3]=$coef[4]=0;
$prev_s = compute_z();
my $min_s = $prev_s;
my @min_coef;
my $min_zscore;
my $min_gap;
my $min_gap1;
my $min_gap10;
my $min_eref;

my $min_T = 0.00001;
my $KT = 0.2;
printf("Init: %f %f : %f %f %f %f e_ref: %f score: %f gaps: %f %f %f\n", $KT, $min_s, $coef[1],$coef[2],$coef[3],$coef[4], $e_ref, $zscore, $gap, $gap1, $gap10);
for(my $nt=0; $nt<200; $nt++){
  $KT = $min_T + ($KT-$min_T)*0.9;
  $index=0;
  while($index<1000){
    my $i = rand(4)+1;
    my $c = rand($max_c*2)-$max_c;
    my $prev_coef = $coef[$i];
    $coef[$i] = $c;
    my $new_s = compute_z();
    
    my $ds = $new_s-$prev_s;
    my $prob = exp(-$ds/$KT);
    
    if($new_s <= $prev_s || rand() < $prob){
      #printf("%d %f: %f %f %f %f\n",$index, $new_s, $coef[1], $coef[2], $coef[3], $coef[4]);
      $prev_s = $new_s;
    }
    else{
      $coef[$i] = $prev_coef;
    }
    if($prev_s < $min_s){
      $min_s = $prev_s;
      for(my $k=0; $k<5; $k++){
	$min_coef[$k] = $coef[$k];
      }
      $min_zscore = $zscore;
      $min_gap = $gap;
      $min_gap1 = $gap1;
      $min_gap10 = $gap10;
      $min_eref = $e_ref;
    }
    $index++;
  }
  $prev_s = $min_s;
  for(my $k=0; $k<5; $k++){
    $coef[$k] = $min_coef[$k];
  }
  
  printf("%d %f %f : %f %f %f %f e_ref: %f score: %f gaps: %f %lf %lf\n", $nt, $KT, $min_s, $min_coef[1],$min_coef[2],$min_coef[3],$min_coef[4], $min_eref, $min_zscore, $min_gap, $min_gap1, $min_gap10);
}



sub compute_z(){
  my $p_ave=0;
  my $p_ave2=0;
  my $max = -999999999;

  my @array1;
  for(my $i=0; $i<$np; $i++){
    my $tmp=0;
    for(my $k=0; $k<5; $k++){
      $tmp+= $positive[$i][$k]*$coef[$k];
    }
    if($tmp > $max){
      $max = $tmp;
    }
    $array1[$i] = $tmp;
    $p_ave += $tmp;
    $p_ave2 += $tmp*$tmp;
  }
  $p_ave /= $np;
  $p_ave2 /= $np;
  $p_ave2 = sqrt($p_ave2 - $p_ave*$p_ave);
  
  
  
  my $n_ave=0;
  my $n_ave2=0;
  my $min = 999999999;
  my @array2;
  for(my $i=0; $i<$nn; $i++){
    my $tmp=0;
    for(my $k=0; $k<5; $k++){
      $tmp+= $negative[$i][$k]*$coef[$k];
    }
    if($tmp < $min){
      $min = $tmp;
    }
    $array2[$i] = $tmp;
    $n_ave += $tmp;
    $n_ave2 += $tmp*$tmp;
  }
  $n_ave /= $nn;
  $n_ave2 /= $nn;
  $n_ave2 = sqrt($n_ave2 - $n_ave*$n_ave); 
  
  my @sorted_1 = sort { $a <=> $b } @array1;
  my @sorted_2 = sort { $a <=> $b } @array2;
  
  #printf("%f %f; %f %f\n", $sorted_1[0], $sorted_1[$np-1], $sorted_2[0], $sorted_2[$nn-1]);
  my $p_ave10=0;
  my $cut = 0.1;
  for(my $i=0; $i<$np*$cut; $i++){
    $p_ave10 += $sorted_1[$np-1-$i];
  }
  $p_ave10 /= int($np*$cut);
  my $n_ave10=0;
  for(my $i=0; $i<$nn*$cut; $i++){
    $n_ave10 += $sorted_2[$i];
  }
  $n_ave10 /= int($nn*$cut);
  
  $c1=0;
  $e_ref = 0;
  for(my $i=0; $i<$np; $i++){
    my $a = int($nn*($i+1)/$np);
    if($sorted_2[$a] > $sorted_1[$np-1-$i]){
      $c1 = $i/$np;
      $e_ref = $sorted_2[$a];
      last;
    }
  }
  
  $zscore = ($p_ave - $n_ave)/($p_ave2+$n_ave2);
  $gap1 = ($max-$min)/($n_ave-$p_ave);
  $gap10 = ($p_ave10 - $n_ave10)/($n_ave-$p_ave);
  $gap = $c1;
  #printf("%f %f %f %f\n", $p_ave, $p_ave2, $n_ave, $n_ave2);
  return $zscore + $gap10 + $gap1*0.1 + $gap*2;
}
