#!/usr/bin/perl

die "Usage: $ARGV[0] <basename> <rho_start> <rho stop> <rho step>$#ARGV\n" 
  if $#ARGV < 3;

$basename=$ARGV[0];
$start=$ARGV[1];
$stop=$ARGV[2];
$rhostep=$ARGV[3];
$script=$basename . ".as";
$script2=$basename . "2.as";
$likedata=$basename . "_like.dat";

open LIKEDATA, ">$likedata" or die "Can't open $likedata\n";
print LIKEDATA  "Rho LogLike RhoAdj Deviance\n";

for($rho=$start;$rho<=$stop;$rho+=$rhostep){
  open SCRIPT , "< $script" or die "Can't open $script\n";
  open OUTSCRIPT,"> $script2" or die "Can't open $script2\n";

  while(<SCRIPT>){
    s/RHOVAR/$rho/;
    print  OUTSCRIPT $_;
  }

  system("ASReml -n $script2 > log");

  open RESULTS, "<$script2" . "r" or die "Can't open $script2" . "r/n";
  while(<RESULTS>){
    chomp();
    @parts=split/ +/;
    $rhoadj=$parts[6] if $parts[2]=~/RhoAdj/;
    $dev=$parts[6] if $_=~/Deviance from GLM fit/;
  }

  $loglike=$rhoadj-$dev/2.;
  if($rho==$start | $loglike> $maxlike){
    $maxlike=$loglike ;
    $maxrho=$rho;
  }
  print LIKEDATA "$rho $loglike $rhoadj $dev\n";

}
print "MLE rho=$maxrho with logLike=$maxlike\n";