#!/usr/bin/perl -w
use strict;
my $dir      = "Rfam";
my $outdir   = "./results";
my $seedfile = "seeds.txt";
my $l        = 32;
goto IV;


I:

# foreach seed
open S,$seedfile;
while (<S>) {

    #---------------------#
    # I) Compute features #
    #---------------------#
    chomp $_;
    if ($_ eq ""){
	next;
    }
    my $seeds     = "$_";
       $seeds =~ /^\Q${seeds}\E/;
    my $features    = $outdir."/features.".${seeds}.".Rfam.ft";
    my $train       = $outdir."/train.".${seeds}.".Rfam.ft";
    my $test        = $outdir."/test.".${seeds}.".Rfam.ft";
    my $model       = $outdir."/model.".${seeds}.".Rfam.svm";
    my $predictions = $outdir."/predictions.".${seeds}.".Rfam.res";
    
    # compute the features for each fasta sequence of the learning set
    opendir(DIR,$dir);
    my @files = readdir(DIR);
    closedir(DIR);
    system("rm -f \"${features}\"");

    # foreach fasta file
    foreach my $file (@files) {
	if ($file eq "." || $file eq ".."){
	    next;
	}
	
	# get and print the class header
	open F, ">> ${features}" or die "Can't open ${features} : $!";    
	my $header = "";
	if ($file =~ m/.*RF[0]*([1-9][0-9]*)_.*$/) {
	    $header = $1; 
	} else {
	    $header = $file;
	}
	print F $header;
	close F;
	
	# call the C++ program that computes features
        system("./dna_features \\\"${seeds}\\\" \"${dir}/${file}\" >> \"${features}\" 2> /dev/null\n");
    }

    #------------------------#
    # II) Learn and Classify #
    #------------------------#

    system("./split  \"${features}\" \"${train}\" \"${test}\"\n");
    system("./svm_multiclass_learn -c 0.1 \"${train}\" \"${model}\"\n");
    system("./svm_multiclass_classify     \"${test}\"  \"${model}\"  \"${predictions}\" 1> \"${predictions}.stdout\" 2> \"${predictions}.stderr\"\n");
    # added to avoid too much space usage (must be removed you want to keep the classifier models)
    system("rm -f \"${features}\" \"${train}\" \"${test}\" \"${model}\" \"${predictions}\"\n");
}

III:

#------------------------#
# III) Print results     #
#------------------------#

system("rm -f results.dat");

# foreach seed
open S,$seedfile;
while (<S>) {
    chomp $_;
    if ($_ eq ""){
	next;
    }
    my $seeds       = "$_";
    my $features    = $outdir."/predictions.".${seeds}.".Rfam.res.stdout";
    # display seeds and error rate of svm
    system("cat ${features} | grep \"Zero.one-error on test set:\" | sed 's/Zero.one-error on test set: \\(\[0-9\.\]*\\)%.*\$/\t\\1/' | tr \"\\n#-\" \"\\t1*\" >> results.dat");
    #sleep(1);
    # and show the multihit sensitivity
    for (my $x = 1; $x < $l-7; $x++) {
      system("./iedera -spaced -l ${l} -y $x -f 0.3,0.7 -m \"${seeds}\" | sed 's/${seeds}/\"${seeds}\"/g' | tr \"\\n#-\" \"\\t1*\" >> results.dat");
      #sleep(1);
    }
    # and show the coverage sensitivity
    for (my $x = 1; $x < $l-7; $x++) {
      system("./iedera -spaced -l ${l} -g $x -f 0.3,0.7 -m \"${seeds}\" | sed 's/${seeds}/\"${seeds}\"/g' | tr \"\\n#-\" \"\\t1*\" >> results.dat");
      #sleep(1);
    }
    # and show the theoretical sensitivity
    system("./iedera -spaced -l ${l}       -f 0.3,0.7 -m \"${seeds}\" | sed 's/${seeds}/\"${seeds}\"/g'  | tr \"#-\" \"1*\" >> results.dat");
}

IV:

system("rm -f results-n1.dat  results-n2.dat results-w3-n1.dat results-w3-n2.dat results-w4-n1.dat results-w4-n2.dat");
system("cat results.dat | grep \",\" > results-n2.dat");
system("cat results.dat | grep -v \",\" > results-n1.dat");
system("cat results-n2.dat | grep -Ev \"1[\\*]*1[\\*]*1[\\*]*1\\\"\" > results-w3-n2.dat");
system("cat results-n2.dat | grep -E  \"1[\\*]*1[\\*]*1[\\*]*1\\\"\" > results-w4-n2.dat");
system("cat results-n1.dat | grep -Ev \"1[\\*]*1[\\*]*1[\\*]*1\\\"\" > results-w3-n1.dat");
system("cat results-n1.dat | grep -E  \"1[\\*]*1[\\*]*1[\\*]*1\\\"\" > results-w4-n1.dat");

#------------------------#
# IV) Print correlation  #
#------------------------#
foreach my $wn ("w3-n1","w3-n2","w4-n1","w4-n2") {
  system("rm -f  main-${wn}.dat");
  system("echo \"# multi-sens(1..".($l-7).")\\tcov-sens(1..".($l-7).")\\tsingle-sens\" >> main-${wn}.dat");
  my $multicol = 0;
  my $covcol   = ($l-8)*2;
  my $singlecol= ($l-8)*4+2; 
  for (my $x = 1; $x < $l-7; $x++) {
    $multicol += 2;
    system("./corr.pl results-${wn}.dat 1 ${multicol}  >> main-${wn}.dat");
    $covcol   += 2;
    system("./corr.pl results-${wn}.dat 1 ${covcol}    >> main-${wn}.dat");
    system("./corr.pl results-${wn}.dat 1 ${singlecol} >> main-${wn}.dat");
    system("echo \"\" >> main-${wn}.dat");
  }      
  system("gnuplot main-${wn}.gp");
  system("gnuplot results-${wn}.gp");
}



