#! /bin/bash
# Last edited on 2008-06-06 23:25:33 by stolfi

export LC_ALL=C
cmd=${0##*/}
usage=( \
  "${cmd} \\" \
  "\n  {WIN_SIZE} \\" \
  "\n  [ -range {BASINI} {BASFIN} ] \\" \
  "\n  < {INFILE}.est \\" \
  "\n  > {OUTFILE}.eps" \
)

# Reads the output of the "pre" classifier (".est" file) Plots the 'K'
# and 'N' probabilities for the bases in the given index range,
# together with the official classification.

basini=0
basfin=999999999
if [ $# -lt 1 ]; then
  error "missing argument" ;
  echo "usage: ${usage[@]}" 1>&2 ; exit 1
fi
winSize="$1"; shift;
while [ $# -gt 0 ]; do
  case "$1" in
    -range ) 
      basini="$2"; basfin="$3"; 
      shift; shift; shift ;;
    -* )
      error "unrecognized option $1" ;
      echo "usage: ${usage[@]}" 1>&2 ; exit 1 ;;
    * ) break;;
  esac;
done

tmp="/tmp/$$"

pltfile="${tmp}.plt"
epsfile="${tmp}.eps"

# Convert the ".est" file to plottable format:

cat \
  | gawk \
      -v ini="${basini}" -v fin="${basfin}" \
      ' BEGIN { n = -1; } \
        /^[\#>]/{ next; } \
        /^[ ]*$/{ next; } \
        (NF==14){ \
          n++; \
          if ((n < ini+0) || (n > fin+0)) { next; } \
          prK = $8; prN = $9; cod = $10; \
          printf "%s %s\n", prK, prN > "/dev/stderr"; \
          gsub(/^[A-Z?][:]/, "", cod); \
          gsub(/[*?\·]$/, "", cod); \
          if (cod !~ /^[KN?]$/) \
            { printf "%d: bad off code = %s\n", FNR, cod > "/dev/stderr"; } \
          prO = "???"; \
          if (cod == "K") { prO = "1.00"; } \
          if (cod == "N") { prO = "0.00"; } \
          printf "%6.4f %6.4f %s\n", prK, prN, prO; \
          next; \
        } \
        // { printf "%d: bad NF = %d\n", FNR, NF > "/dev/stderr"; } \
      ' \
  > ${pltfile}
  
# Plot it:

gnuplot <<EOF
# set encoding iso_8859_1
set title "window size ${winSize}"
set ylabel "Pr(K)"
set xlabel "Base index"
set term postscript eps color "TimesRoman" 14
set output "${epsfile}"
set size 2,0.4
set yrange [-0.05:+1.05]
plot \
  "${pltfile}" using 0:1 title "K" with histeps lt 1, \
  "${pltfile}" using 0:3 title "N" with histeps lt 7
quit
EOF

ghostview ${epsfile}
cat ${epsfile}
/bin/rm -f ${pltfile} ${epsfile}