#### damatfrc.m
#### Francesco Potortì 2007

## Unchanged parameters
Tp = 252;			# propagation delay [ms]
pktsize = 1316*8/1000;		# size of a packet [kbit]
mindama = .044;			# minimum DAMA rate is 44 kb/s [Mb/s]
maxdama = 144*mindama;		# maximum DAMA rate is 6.3 Mb/s [Mb/s]
Tf_btp = 3;			# allocation periods per BTP
plotres = 500;			# number of points to plot

## Parameters read from stdin
#Tq = 450;			# QSD delay [ms]
#offset = 820;			# when first BTP is received [ms]
#RTTm = 777;			# mean RTT is 2*propagation+frame [ms]
#full = 256;			# streaming rate [kb/s]
#slow_increase = true;		# TFRC uses slow increase
#decay = 2;			# constant for decreasing DAMA allocation [frames]
#Tb = 820;			# superframe duration [ms]
#Ta_btp = 2;			# allocation delay [superframes]

## Read parameters from stdin
if (length(argv)==0 || strcmp(argv(){1}, "cgi")==0)
  error("First argument should be \"cgi\""); endif
disp(["Content-Type: text/html\nConnection: close\n\n",
      "<!DOCTYPE HTML PUBLIC \"-//W3C//DTD HTML 4.01 Transitional//EN\"",
      " \"http://www.w3.org/TR/html4/loose.dtd\"><HTML><HEAD>",
      "<TITLE>Damatfrc.m Output</TITLE>",
      "<META HTTP-EQUIV=\"Content-Type\" CONTENT=\"text/html; charset=iso-8859-1\">",
      "</HEAD><BODY><DIV ALIGN=center><TABLE CELLPADDING=15>\n",
      "<TR><TD VALIGN=top><P>Input parameters<P><PRE>\n"]);
params = char(strsplit(fgetl(stdin),"&"));
for i=1:rows(params)
  param = deblank(params(i,:));
  if (index(param, "("))
    printf("Illegal CGI parameter \"%s\"\n", param);
    error("Illegal CGI parameter \"%s\"", param); endif
  eval(param);
endfor
server = getenv("SERVER_NAME");
referer = getenv("HTTP_REFERER");
hoststart = findstr(referer, server);
dirstart = hoststart + index(referer(hoststart:end), "/") - 1;
dirend = rindex(referer, "/");
dirname = referer(dirstart:dirend);
switch imgfmt
  case 2 fmt = "svg";
    imgtype = "svg+xml";
  otherwise fmt = "png";
    imgtype = "png"
endswitch
imgh = floor(max(48, min(48000, imgh))); # img height
imgw = floor(1.5*imgh);		# img width
imgname = sprintf("%sdamatfrc-%d.%s", dirname, getpid(), fmt);
disp("</PRE><TD VALIGN=top><P>Trace of QS feedback reports<P><PRE>");

## Parameter initialisation
Ta = Ta_btp*Tb;			# allocation delay
Tf = Tb/Tf_btp;			# allocation period
Tr = Tb;			# measurement interval
RTTm = floor(max(2*Tp, min(2*Tp+Tf, RTTm))); # mean RTT
Tq = floor(max(0, min(Tr+Ta-RTTm, Tq))); # QSD delay
offset = floor(max(2, min(Tb, offset))); # time of first BTP ending
full = max(1, min(1000, full))/1000; # streaming rate [Mb/s]
alpha = 1-1/max(1, min(100, decay)); # DAMA deallocation decay constant

qsreq = max(.08,.01*2^ceil(log(full/.01)/log(2))); # QS request rate [Mb/s]
Da1 = offset+Ta;		# time the first DAMA allocation arrives
## The three values Tqopt, thmaxt and thmaxd are approximations valid
## in the case when the delay grows initially and then begins to go
## down. The minimum Tq values for this to happen is called the
## optimum case.
Tqopt = max(1, ceil((2*(Tr+Ta)-3*RTTm)/5)); # optimum Tq (QSD delay) value
mintfrc = pktsize/RTTm;	# minimum TFRC allowed rate [Mb/s]
Frstep = RTTm+Tq;		# Feedback report step
Fr1 = Frstep;	  # time the first feedback report (QS approval) arrives
thmaxt = ((Tr+Ta)-(RTTm+Tqopt))*mindama/full; # when most delayed packet is sent
thmaxd = (Tr+Ta)+RTTm/2-thmaxt; # delay of most delayed packet

tit1 = sprintf("Tb=%d, RTTm=%d[%d,%d], offset=%d[2,%d], alpha=%g, sl_inc=%d, Tq=%d",
	       Tb, RTTm, 2*Tp, 2*Tp+floor(2*Tf), offset, Tr, alpha, slow_increase, Tq);
tit2 = sprintf("Tqopt=%g, th maxdelay=%g", Tqopt, thmaxd);

#############################################################################
########################## Now compute ######################################

## compute cumulative output given previous output, rate,
## upstream output including the previous
function out = outfromrate (prevout, rate, upoutp)
  global err;
  if (abs(prevout)/upoutp(1) > 1+err)
    error("prevout:%g > upoutp(1):%g", prevout, upoutp(1))
  endif
  upout = upoutp(2:end);
  uprate = diff(upoutp);
  lowrate = find(rate./uprate < 1+err);
  maxloop = 10; for loop = 1:maxloop
    out = prevout+cumsum(rate);
    f = find(out./upout > 1+err);
    if (isempty(f))
      return;
    %% else printf("f=[%d..(%d)..%d]\n", f(1), length(f), f(end));
    endif
    l = lowrate(lowrate/f(1) > 1+err);
    if (isempty(l))
      intv = f(1):length(rate);
    else
      %% printf("l=[%d..(%d)..%d]\n", l(1), length(l), l(end));
      intv = f(1):l(1)-1;
    endif
    %% disp([intv(1), intv(end)]);
    rate(intv) = uprate(intv);
  endfor
  error("looping in outfromrate more than %d times", maxloop);
endfunction

## time
t = [0:99999]';		# [ms]
global err = eps*10*length(t);

## output from stream generator
streamrate = full * ones(size(t));
stream = full * t;

## max allowed TFRC rate
tfrcrate = -ones(size(t));	# max rate allowed by TFRC
Frintv = 1:Fr1-1;		# from 0 to the first feedback report
tfrcrate(Frintv) = pktsize/Fr1;
Frintv = Fr1:Fr1+2*Frstep-1;	# from the 1st to the 3rd feedback rep.
tfreport = qsreq;		# QS approval
printf("tfreport: %g received at %d (QS approval), next one at %d (ignored)\n",
       tfreport, Fr1, Fr1+Frstep);
tfrcrate(Frintv) = tfreport;	# not doubled because it is QS
oldrate = newrate = tfreport;
Frintv = 1:Frintv(end);	# interval to compute tfrc on
%% printf("Frintv = %d:%d\n", Frintv(1), Frintv(end));

## output from TFRC
tfrc = -ones(size(t));	# output from TFRC
tfrc(Frintv) = outfromrate(-tfrcrate(1), tfrcrate(Frintv),
			   [-streamrate(1); stream(Frintv)]);
Frintv = Frintv(end)+(1:Frstep); # next interval to compute
%% printf("Frintv = %d:%d\n", Frintv(1), Frintv(end));
## the 3rd feedback report is based on the traffic arriving at the
## receiver in the interval Fr1+Frstep*(1:2)-RTTm/2, or, equivalently,
## on the traffic exiting from DAMA one RTTm/2 before.
qst = Fr1+Frstep-RTTm;	# start of interval
qen = qst+Frstep-1;		# end of interval

## DAMA allocation
damarate = -ones(size(t));	# DAMA allocation
Daintv = 1:Da1-1;		# from 0 to the first allocation
damarate(Daintv) = mindama;
Daintv = Da1:Da1+Tb-1;		# from the 1st to the 2nd allocation
dst = 1; den = offset-1;	# interval the 2nd allocation is based on
if (den>0)
  if (tfrc(den)<0) error("tfrc invalid in %d", den); endif
  traffic = (tfrc(den)-tfrc(dst))/Tr;
else traffic = 0;
endif
allocation = max(mindama, qsreq + traffic);
%% printf("allocation: %g at %d, from tfrc out in %d:%d\n", allocation, Daintv(1), dst, den);
damarate(Daintv) = min(maxdama, allocation);
Daintv = 1:Daintv(end);	# interval to compute dama on
%% printf("Daintv = %d:%d\n", Daintv(1), Daintv(end));

## output from DAMA
dama = -ones(size(t));		# output from DAMA
dama(Daintv) = outfromrate(-damarate(1), damarate(Daintv),
			   [-tfrcrate(1); tfrc(Daintv)]);
Daintv = Daintv(end)+(1:Tb);	# next interval to compute
%% printf("Daintv = %d:%d\n", Daintv(1), Daintv(end));
dst = den+1; den += Tb; # interval the next DAMA allocation is based on

## Alternately compute tfrcrate and damarate until the end of time.
lastt = length(t)-4*Ta;
do
  if (tfrc(qen)!=-1 && damarate(qen)!=-1)
    dama(qst:qen) = outfromrate(dama(qst-1), damarate(qst:qen), tfrc(qst-1:qen));
    if (stream(qen)-dama(qen)<pktsize/50) lastt = qen; endif
    tfreport = (dama(qen)-dama(qst))/(qen-qst); # feedback report
    oldrate = newrate;
    if (oldrate>2*tfreport)	# report requests reduction
      oldrate = newrate = max(mintfrc, tfreport);
    else
      ## Should compute the current rate, like this, but we don't have tfrc yet:
      ## currate = tfrc(Frintv(1)-1) - tfrc(Frintv(1)-2);
      ## however, in our case approximating with oldrate is exact.
      newrate = min(2*oldrate, max(mintfrc, 2*tfreport));
    endif
    tfrcrate(Frintv) = newrate;
    printf("tfreport: %g at %d, from dama out in %d:%d, newrate: %g\n",
	   tfreport, Frintv(1), qst, qen, newrate);
    if (slow_increase && newrate>oldrate)
      tfrcrate(Frintv) += (newrate-oldrate)*(Frintv'-Frintv(end)-1)/(range(Frintv)+1);
    endif
    Frintv += Frstep;		# next interval to compute
    qst += Frstep; qen += Frstep;	# next interval the feedback report is based on
    %% printf("Frintv = %d:%d\n", Frintv(1), Frintv(end));
  elseif (tfrcrate(den)!=-1)
    tfrc(dst:den) = outfromrate(tfrc(dst-1), tfrcrate(dst:den), stream(dst-1:den));
    prevall = allocation;
    traffic = (tfrc(den)-tfrc(dst))/(den-dst);
    allocation = mindama*max(1, ceil(traffic/mindama));
    %% printf("allocation: %g at %d, from tfrc out in %d:%d\n",
    %% allocation, Daintv(1), dst, den);
    if (allocation<prevall) allocation += (prevall-allocation)*alpha; endif
    damarate(Daintv) = min(maxdama, allocation);
    Daintv += Tb;		# next interval to compute
    dst += Tb; den += Tb;	# next interval the DAMA allocation is based on
    %% printf("Daintv = %d:%d\n", Daintv(1), Daintv(end));
  elseif (tfrcrate(qen)!=-1)
    tfrc(dst:qen) = outfromrate(tfrc(dst-1), tfrcrate(dst:qen), stream(dst-1:qen));
  else
    error("cannot compute next interval");
  endif
until (qen>lastt || Frintv(end)>t(end) || Daintv(end)>t(end))
disp("</PRE></TABLE>");

## Computation finished; prepare plotting data
plotstep = ceil(lastt/plotres);
dest = prepad(dama(1:end-floor(RTTm/2)), length(t)); # seq. no. at destination
intv = (dest>0 & t<lastt & mod(t,plotstep)==0);
delayx = dest(intv)/full;
delayy = t(intv) - delayx;
x = t(1:plotstep:lastt);
y = [stream, tfrc, dama, dest](1:plotstep:lastt,:);
Frtime = Fr1:Frstep:lastt;	# times feedback reports are received
Datime = Da1:Tb:lastt;	# times DAMA allocations are received
titd = sprintf("maxdelay=%g", max(delayy));
root = getenv("DOCUMENT_ROOT");
imgfilename = strcat(root, imgname);
unlink(imgfilename);

## Do the plot
graphics_toolkit gnuplot
#gnuplot_binary("/usr/local/bin/gnuplot4.2-nox");
set (0, "defaultaxesfontname", "Helvetica");	# bug #710272
titstr = strcat("One-way delay and cumulative traffic",
		" at the output of: generator, TFRC, DAMA, receiver\n",
		tit1, "\n", tit2, ", ", titd);
figure ("visible", "off");
plot(x, y(:,1), "-b;;", x, y(:,2), "-g;;", x, y(:,3), "-c;;", x, y(:,4), "-m;;",
     Frtime, tfrc(Frtime), "*g;feedback reports;",
     Datime, dama(Datime), "ob;DAMA allocations;",
     delayx, delayy, "-r;delay;",
     "markersize", 1);
title(titstr);
xlabel("Time from start of generator [ms]");
ylabel("Delay [ms] and cumulative traffics [kbit]");
legend('location','northwest');
print(imgfilename, sprintf("-d%s", fmt), sprintf("-S%d,%d", imgw, imgh));

## Output the final HTML and wait for gnuplot to finish
filesize = 0;
for i=1:20			# 10s timeout on gnuplot
  [s err] = stat(imgfilename);
  if (err==0)
    oldfilesize = filesize;
    filesize = s.size;
    if (filesize>0 && filesize==oldfilesize)
      close;
      pause(0.5);		# so that 1s elapses from img ctime to atime
      printf("<OBJECT TYPE=\"image/%s\" DATA=\"%s\" WIDTH=%d HEIGHT=%d></OBJECT>\n",
	     imgtype, imgname, imgw, imgh);
      disp("</DIV></BODY></HTML>");
      exit;			# finished with CGI
    endif
  endif
  # printf("<p>Waiting for Gnuplot...\n");
  pause(0.5);
endfor
printf("<P>Gnuplot timeout</BODY></HTML>\n");
error("Gnuplot timeout");
