#!/bin/csh

setenv _POSIX2_VERSION 199209

# Copyright (C) 2005 Andrew E Firth, University of Otago, Dunedin, 
# New Zealand, aef(at)sanger.otago.ac.nz
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License (version 2) as
# published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
# 02110-1301, USA.

#Redraws MLEnuc plot.

limit stacksize unlimited

@ argc = `echo $argv | awk '{print NF}'`
if ($argc != 8) then
  echo "Usage: redoMLEnuc newdir halfwindow base1 base2 grids2 count threshold newurl<br>"
  exit 1
endif

set newdir = $1
cd $newdir
source mlogd.param
source MLEnuc.params

if (! $download) then
  setenv PATH ${PATH}:/home/aef/bin
else
  setenv PATH ${PATH}:${HOME}/bin
endif

set halfwindow = $2
@ base1 = $3
@ base2 = $4
set grids = $5
set count = $6
set threshold = $7
set newurl = $8

if ($base1 > $base2) then
  @ temp = $base1
  @ base1 = $base2
  @ base2 = $temp
endif

set lambdasum2 = `head -2 tree.info | tail -1 | awk '{print $1*0.25}'`
set threshold = `echo $threshold | awk '{if (0.+$1<0.01) {print "0.01"} else {print $1}}' | awk '{if (0.+$1>0.99) {print "0.99"} else {print $1}}'`
set threshold2 = `echo $threshold $lambdasum2 | awk '{print $1*$2}'`

if (1 == $grids) then
  set tick = ""
else
  set tick = "#"
endif

echo -n "Preparing plot data...<br>"

#Make running mean tracks for ref-nonref pairs.
awk '{print $2}' ref.MLElog | sort | uniq > temp1
rm -f ref.runmean
touch ref.runmean
foreach nrseq (`cat temp1`)
  awk '{if ($2=="'"$nrseq"'") print $0}' ref.MLElog > temp2
  runmean2 temp2 $halfwindow $circular >> ref.runmean || exit 1
end

touch ref.runmean
rm -f temp[12]

#Make running mean tracks for regions of treesum with summed divergence
#  >= $threshold * $lambdasum2.
awk '{if (0.+$3>=0.+"'"$threshold2"'"&&0.+$3>0.) print $1,$2,0,$4*"'"$lambdasum2"'"/$3,$5*"'"$lambdasum2"'"/$3}' tree.treesum > temp1
runmean2 temp1 $halfwindow $circular > tree.runmean || exit 1
touch tree.runmean
rm -f temp[12]

echo "...Done<br><br>"

echo -n "Making plots...<br>"

#Process mlogd.R and plot.
set date = `date +%y-%m-%d`
set window = `echo $halfwindow | awk '{print 1+$1*2}'`
awk '{printf "%14.6f\n",sqrt(($5-$4)*($5-$4))}' ref.runmean > temp1
set yyyscale = `minmax temp1 | awk '{if (0.+$2>0.) {print 0.5/$2} else {print "0"}}'`
rm -f temp1
awk '{printf "%14.6f\n",$5-$4}' tree.runmean > temp1
echo 0 >> temp1
set y2min = `minmax temp1 | awk '{print $1}'`
set y2max = `minmax temp1 | awk '{print $2}'`
rm -f temp1

#Check for empty files to edit out of R script to stop R from barfing.
@ w = `wc -l ref.runmean | awk '{print $1}'`
if ($w == 0) then
  set kkk6 = "#"
else
  set kkk6 = ""
endif
@ w = `wc -l tree.runmean | awk '{print $1}'`
if ($w == 0) then
  set kkk7 = "#"
else
  set kkk7 = ""
endif

cat ../../SCRIPTS/mlogd.R | \
  sed 's/MLEnuc/MLEnuc.'$count'/g' | \
  sed 's/MLEscatter/MLEscatter.'$count'/g' | \
  sed 's/MLEtrack/MLEtrack.'$count'/g' | \
  sed 's/nnn/'$nseqs'/g' | sed 's/ddd/'$date'/g' | \
  sed 's/bbb1/'$base1'/g' | sed 's/bbb2/'$base2'/g' | \
  sed 's/bbbq1/'$base1'/g' | sed 's/bbbq2/'$base2'/g' | \
  sed 's/www/'$window'/g' | sed 's/lll/'$lambdasum'/g' | \
  sed 's/orf1x1/'$orf1x1'/g' | sed 's/orf1x2/'$orf1x2'/g' | \
  sed 's/orf2x1/'$orf2x1'/g' | sed 's/orf2x2/'$orf2x2'/g' | \
  sed 's/yyyscale/'$yyyscale'/g' | sed 's/fff/'$ntmlesum'/g' | \
  sed 's/#vvv0/'$tick'/g' | \
  sed 's/#vvv1/'$vvv1'/g' | sed 's/#kkk0/'$kkk0'/g' | \
  sed 's/#kkk1/'$kkk1'/g' | sed 's/#kkk2/'$kkk2'/g' | \
  sed 's/#kkk3/'$kkk3'/g' | sed 's/#kkk4/'$kkk4'/g' | \
  sed 's/#kkk5/'$kkk5'/g' | sed 's/#kkk6/'$kkk6'/g' | \
  sed 's/#kkk7/'$kkk7'/g' | sed 's/#kkk9//g' | \
  sed 's/yyy2max/'$y2max'/g' | sed 's/yyy2min/'$y2min'/g' | \
  sed 's/sss3/'$sss3'/g' | sed 's/sss1/'$sss1'/g' | \
  sed 's/sss2/'$sss2'/g' | sed 's/@/"/g' | \
  sed 's/%sss%/'$refseq'/g' | sed 's/%ggg%/'$title'/g' > plot1.R
R --no-save --slave < plot1.R >>& errorlog.txt
epstopdf MLEnuc.$count.eps 
epstopdf MLEscatter.$count.eps 

echo "...Done<br><br>"

chmod a+rw ./*
#Write rest of webpage.
set id9 = `echo $newurl | sed 's/.*\///'`
echo "<hr><b>Nucleotide-by-nucleotide plot:</b> (download "
echo "<a href=$newurl/MLEnuc.$count.png>png</a>, "
echo "<a href=$newurl/MLEnuc.$count.eps>eps</a>, "
echo "<a href=$newurl/MLEnuc.$count.pdf>pdf</a>).<br>"
echo "<b>Zoomed-in plot:</b> (download "
echo "<a href=$newurl/MLEtrack.$count.png>png</a>).<br>"
if (! $download) then
  echo "<a href=redo_MLEnuc.pl?id=$id9&count=$count>Redraw</a>"
  echo "these plot again.  (Don't just use the back button in your browser, "
  echo "or the plot may not refresh.)<br>"
  echo "Return to <a href=mlogd.pl>MLOGD server</a>.<br>"
endif
echo "<img src="'"'"$newurl/MLEnuc.$count.png"'"'" "
echo "alt="'"'"nuc-by-nuc MLE plot (png)"'"'"><br><br>"
echo "<img src="'"'"$newurl/MLEtrack.$count.png"'"'" "
echo "alt="'"'"nuc-by-nuc MLE plot (png) - too long to plot - try Save Link Target As"'"'"><hr>"

chmod a+rw ./*