#! /bin/bash
# Last edited on 2014-03-18 20:19:00 by stolfilocal

# Scatterplots of slumber volume ratio and daily volume.
# 
# Arguments: {TITLE} {PMIN} {PMAX} {WH0} {WH1} {WH2} {REFS} {FNAME}
# 
# Where 
#   
#   {TITLE} is the plot's title;
#   {WH0,WH1,WH2} are coefficients to mix {Vh0,Vh1,Vh2} into the slumber volume {Vh}.
#   {REFS} is the soft threshold for the {Vh/Vd} ratio.
#   {FNAME} is a data file with one line per slumber time in format
#      "{DATE} {TIME} | {NDAY} | {QL} | {HUBI.VH0} | {HUBI.VH1} | {HUBI.VH2} | {HUBI.VD}"

title="$1"; shift  # Plot title.
wh0="$1"; shift    # Scaling for {Vh0}.
wh1="$1"; shift    # Scaling for {Vh1}.
wh2="$1"; shift    # Scaling for {Vh2}.
refS="$1"; shift   # Soft {Vd/Vd} threshold.
fname="$1"; shift  # Input file name. 

show="SHOW"

tmp="/tmp/$$"

fullimg="${tmp}-full.png"

color=( '#0022ff' '#ff0000' '#008800' '#8800dd' '#dd4400' '#0066ff' )

export GDFONTPATH=.:..

gnuplot <<EOF
set term png size 1600,1600 font "courbd,20"
set output "${fullimg}"

# ----------------------------------------------------------------------
# Function definitions:

# Volume scaling weights:
wh0 = ${wh0} # Scaling for {Vh0}. 
wh1 = ${wh1} # Scaling for {Vh1}. 
wh2 = ${wh2} # Scaling for {Vh2}. 

refS = ${refS} # {Vh/Vd} threshold.

# Mean slumber volume, excluding zero volumes (NaN if all weights*volumes are zero): 
vhmed(Vh0,Vh1,Vh2) = (wh0*Vh0 + wh1*Vh1 + wh2*Vh2)/((Vh0==0 ? 0 : wh0) + (Vh1==0 ? 0 : wh1) + (Vh2==0 ? 0 : wh2))

minwt = 0.0   # Min weight.
minps = 1.0   # Min point size.
refps = 7.0   # Max point size.
refwt = 1.0   # Weight for point size = {refps}.

load "plot_slumber_funcs.gnuplot"

vmplot(ih0,ih1,ih2) = vhmed(column(ih0),column(ih1),column(ih2))
rvmplot(ih0,ih1,ih2,id) = vmplot(ih0,ih1,ih2)/column(id)
vdplot(id) = column(id)

gauss(z) = (z > 80 ? 0.0 : exp(-0.5*(z)**2)) # Gaussian bell (1 at 0).
dcexp(z) = (z > 100 ? 0.0 : exp(-z))         # Decaying exponential (1 at 0).

# Weight of point computed from Vh,Vd:
weight(Vh,Vd) = ((wh0+wh1+wh2 == 0) || (Vh == 0) || (Vd == 0) || (refS == 0) ? 1.0 : gauss((Vh/Vd)/refS))
weight3(Vh0,Vh1,Vh2,Vd) = weight(vhmed(Vh0,Vh1,Vh2),Vd)

wpsQL(iq) = wps(column(iq)) # Point size from QL column number.
wpsVV(ih0,ih1,ih2,id) = wps(weight3(column(ih0),column(ih1),column(ih2),column(id)))

# ----------------------------------------------------------------------
# Axes and stuff:

set timefmt "%Y-%m-%d %H:%M:%S" # For {timecolumn}.
set lmargin 9
set rmargin 4
set tmargin 2
set bmargin 5

set xtics mirror
set grid xtics lt 1 lw 3 lc rgb '#ffddaa', lt 1 lw 1.5 lc rgb '#ffddaa'
set grid mxtics
set xlabel  "${title} - mean daily volume Vd (kBTC)"

set ytics mirror
set grid ytics lt 1 lw 3 lc rgb '#ffddaa', lt 1 lw 1.5 lc rgb '#ffddaa'
set grid mytics
set ylabel  "${title} - slumber volume ratio Vh/Vd (x 1000)"

set xrange [-0.01:]
set yrange [:]

plot \
  (refS) title "split" with lines lt 1 lw 1.5 lc rgb '#dd0077', \
  "< egrep -e '^20' ${fname}" using (vdplot(14)):(rvmplot(8,10,12,14)):(wpsQL(6)) title "hand weight" \
    with points pt 6 ps variable lt 1 lw 2.0 lc rgb '#0077dd', \
  "" using (vdplot(14)):(rvmplot(8,10,12,14)):(wpsVV(8,10,12,14)) title "computed weight" \
    with points pt 2 ps variable lt 1 lw 2.0 lc rgb '#dd0077'
    
quit
EOF

if [[ -s ${fullimg} ]]; then
  outimg="${tmp}.png"
  convert ${fullimg} -resize '50%' ${outimg}

  if [[ "/${show}" == "/SHOW" ]]; then
    display ${outimg}
  fi

  cat ${outimg}
fi

rm -fv ${tmp}{-*,}.*