#!/usr/bin/octave -qf

## This script implements the models and produces the plots shown in the
## paper:
##
## S. Ferretti, V. Ghini, M. Marzolla, F. Panzieri, Modeling the Energy
## Consumption of Multi-NIC Communication Mechanisms, proc. IEEE Online
## Conference on Green Communications (GreenCom), 25--28 september 2012,
## pp. 104—109, ISBN 978-1-4799-0395-5 
## doi http://dx.doi.org/10.1109/GreenCom.2012.6519624
##
## The script is written in GNU Octave (http://www.octave.org/).
## To use this script you need to install the Octave Queueing toolbox
## available here: http://www.moreno.marzolla.name/software/queueing/
##
## This script was written by Moreno Marzolla <moreno.marzolla (at)
## unibo.it> on jan-jul 2012. This script is released by the author
## according to the Creative Commons CC0 1.0 Universal license
## http://creativecommons.org/publicdomain/zero/1.0/

1; # not a function file

figure(1,"visible","off");
figure(2,"visible","off");

more off;

## Global plot parameters
set(1,"paperorientation","landscape");
papersize=[5 4] * 0.85; 
margin=[0 0];
set(1,"papersize",papersize);
set(1,"paperposition", [margin papersize-margin*2]);

set(2,"paperorientation","landscape");
papersize=[5 4] * 0.4; 
margin=[0 0];
set(2,"papersize",papersize);
set(2,"paperposition", [margin papersize-margin*2]);

##############################################################################
## Builds the infinitesimal generator matrix for a single NIC. Parameter
## names are the same used in the paper. States are mapped as follows:
##
## d -> 1
## s -> 2
## c -> 3
## f -> 4
function Q = mkmodel_single( p, aa, bb, gg, mm )
  Q = zeros(4,4);

  Q(1,2) = aa;
  Q(2,1) = (1-p)*bb;
  Q(2,3) = p*bb;
  Q(3,4) = gg;
  Q(4,1) = mm;

  Q -= diag(sum(Q,2));
endfunction

##############################################################################
## Builds the infinitesimal generator matrix for two NICs with SIP.
## Parameter names are the same used in the paper; in this case
## parameters are vectors of two elements. States are mapped as follows:
##
## d  -> 1
## s1 -> 2
## c1 -> 3
## f1 -> 4
## s2 -> 5
## c2 -> 6
## f2 -> 7
function Q = mkmodel_sip( p, aa, bb, gg, mm )
  Q = zeros(7,7);

  Q(1,2) = aa(1);
  Q(2,1) = (1-p(1))*bb(1);
  Q(2,3) = p(1)*bb(1);
  Q(3,4) = gg(1);
  Q(4,1) = mm(1);

  Q(1,5) = aa(2);
  Q(5,1) = (1-p(2))*bb(2);
  Q(5,6) = p(2)*bb(2);
  Q(6,7) = gg(2);
  Q(7,1) = mm(2);

  Q -= diag(sum(Q,2));
endfunction

##############################################################################
function Q = mkmodel_abps( p, aa, bb, gg, mm )

  Q = zeros(16,16);

  Q1 = mkmodel_single(p(1),aa(1),bb(1),gg(1),mm(1));
  Q1(1:5:end) = 0; # set diagonal elements to zero
  Q2 = mkmodel_single(p(2),aa(2),bb(2),gg(2),mm(2));
  Q2(1:5:end) = 0; # set diagonal elements to zero

  for i=1:16
    for j=1:16
      [i1 i2] = ind2sub([4 4], i); # index in Q1
      [j1 j2] = ind2sub([4 4], j); # index in Q2
      if ( i1 == j1 || i2 == j2 ) 
	Q(i,j) = Q1(i1,j1) + Q2(i2,j2);
      endif
    endfor
  endfor  
  if ( mm(1) == 0 )
    Q( sub2ind([4 4],3,4), sub2ind([4 4],3,1) ) = 1;
    Q( sub2ind([4 4],4,3), sub2ind([4 4],1,3) ) = 1;
  endif
  Q -= diag(sum(Q,2));
endfunction

##############################################################################
function ax = makeplot( xrange, yrange, data_single, data_sip, data_abps, title_str )
  figure(1);
  clf;
  hold on;
  surface(xrange, yrange, data_single, "EdgeColor", [1 1 1], "FaceColor", [.2 .2 .2]);
  surface(xrange, yrange, data_sip, "EdgeColor", [1 1 1], "FaceColor", [.4 .4 .4] );
  surface(xrange, yrange, data_abps, "EdgeColor", [1 1 1], "FaceColor", [.6 .6 .6]);
  hold off;
  view(129,26);
  axis("tight");
  ax = axis();
  ax(5) = ax(5)-0.3*(ax(6)-ax(5));
  axis(ax);
  grid("on");
  box("off");  
  title(title_str);
  xlabel("\\gamma_{WiFi}");
  ylabel("\\gamma_{UMTS}");
  ax = axis();
endfunction

##############################################################################
function makeslicedplot(xrange, yrange, 
			axe, # either "x" or "y"
			idx, # index of the slices
			data_single,
			data_sip,
			data_abps, 
			ax # axes 
			)
  hold on;
  for xs=idx
    if ( axe == "x" || axe == "X" )
      XX = xrange(xs)*ones(1,length(xrange));
      YY = yrange;
      ZZ_single = data_single(:,xs);
      ZZ_sip  = data_sip(:,xs);
      ZZ_abps = data_abps(:,xs);
    else
      XX = xrange;
      YY = yrange(xs)*ones(1,length(yrange));
      ZZ_single = data_single(xs,:);
      ZZ_sip  = data_sip(xs,:);
      ZZ_abps = data_abps(xs,:);
    endif
    plot3(XX, YY, ZZ_single, ";UMTS;", "linewidth", 3, \
	  XX, YY, ZZ_sip, ";SIP;", "linewidth", 3, \
	  XX, YY, ZZ_abps, ";ABPS;", "linewidth", 3 );
  endfor
  axis(ax); ## set axes limits
  xl=get(gca,"xlim");
  yl=get(gca,"ylim");
  zl=get(gca,"zlim");
  for xs=idx
    if ( axe == "x" || axe == "X" )
      plot3([1 1 1 1 1]*xrange(xs),[yl(1) yl(1) yl(2) yl(2) yl(1)], [zl(1) zl(2) zl(2) zl(1) zl(1)],"k","linewidth",1)
    else
      plot3([xl(1) xl(1) xl(2) xl(2) xl(1)], [1 1 1 1 1]*yrange(xs), [zl(1) zl(2) zl(2) zl(1) zl(1)],"k","linewidth",1)
    endif
  endfor
  hold off;
  view(129,26);
endfunction

##############################################################################
function makeplot_slices(xrange, yrange, data_single, data_sip, data_abps, ax, name)
  slices = [1 5 10];
  for idx=1:3
    figure(2);
    clf;
    makeslicedplot(xrange, yrange, "y", slices(idx), data_single, data_sip, data_abps, ax);
    if ( idx == 1 )
      legend("show");
    else
      legend("hide");
    endif
    box("off");
    grid("off");
    axis("nolabel");
    print([name "_" num2str(idx) ".eps"], "-deps2", "-tight");
  endfor
endfunction

##############################################################################
##
## MODEL DEFINITION AND EVALUATION
##
##############################################################################

##############################################################################
##
## MODEL PARAMETERS (NIC 1 is UMTS, NIC 2 is wifi)
##

gamma_umts = linspace(0.0004,0.002,10);
gamma_wifi = linspace(0.02,0.16,10);
aa = [1/6.024 1/7.5]; # alpha
bb = [1/1.5 1/1.5]; # beta
gg = [1/550 1/40]; # gamma
mu = [1 1]; # mu
pp = [.99 .9]; # p

T = 200; # Duration of the battery consumption graph

## REWARDS

## Energy consumption rate (W) for state "c" are taken from table 1 of
## Balasubramanian et al. Energy consumption rates for other states are
## not known (at least to me), so they are zero for now.
pwr_max_umts = 12.5/20;
pwr_max_wifi = 7.6/20;

r_energy_umts = pwr_max_umts*[.2 .5 1 .4]; ## FIXME: mie assunzioni
r_energy_wifi = pwr_max_wifi*[.2 .5 1 .4]; ## FIXME: mie assunzioni
r_energy_single = r_energy_umts;
r_energy_sip = \
    [r_energy_umts 0 0 0] + \
    [r_energy_wifi(1) 0 0 0 r_energy_wifi(2:4)];

## Throughput (in Mbps)
bps_umts = 0.2;
bps_wifi = 26;

r_tput_umts = [0 0 bps_umts 0];
r_tput_wifi = [0 0 bps_wifi 0];
r_tput_single = r_tput_umts;
r_tput_sip = [0 0 bps_umts 0 0 bps_wifi 0];

## For ABPS, definition of rewards is a bit tricky. We do that here
r_tput_abps = r_energy_abps = zeros(1,16);
for i=1:4
  for j=1:4
    ii = sub2ind([4 4],i,j);
    if ( i == 3 && j == 3 )
      r_energy_abps(ii) = r_energy_wifi(j) + r_energy_umts(i)*.2; ## FIXME: mie assunzioni
      r_tput_abps(ii) = r_tput_wifi(j);
    else
      r_energy_abps(ii) = r_energy_umts(i) + r_energy_wifi(j);
      r_tput_abps(ii) = r_tput_umts(i) + r_tput_wifi(j);
    endif
  endfor
endfor


##
## END MODEL PARAMETERS
##
##############################################################################

n_umts = length(gamma_umts);
n_wifi = length(gamma_wifi);

## Results
tput_single = tput_sip = tput_abps = \
mttf_single = mttf_sip = mttf_abps = \
energy_single = energy_sip = energy_abps = \
zeros(length(n_umts),length(n_wifi));
batt_single = batt_sip = batt_abps = zeros(1,T);

##
## Dump parameter values to LaTeX-formatted tables which can be included
## in the paper
##
fout = fopen("rewards_power.tex","w");
fprintf(fout,"$\\mathbf{r}_\\text{UMTS}$");
fprintf(fout, " & %.2f", r_energy_umts);
fprintf(fout, "\\\\\n");
fprintf(fout,"$\\mathbf{r}_\\text{Wi-Fi}$");
fprintf(fout, "& %.2f", r_energy_wifi);
fprintf(fout, "\\\\\n");
fclose(fout);

printf("   1...5....10\n");
## Ready to produce the raw data
for i=1:n_umts
  printf("%2d ",i);
  for j=1:n_wifi

    ##
    ## single NIC
    ##

    Q_single = mkmodel_single(pp(1), aa(1), bb(1), gamma_umts(i), mu(1));
    q_single = ctmc(Q_single);
    tput_single(i,j) = dot(q_single,r_tput_single);
    energy_single(i,j) = dot(q_single,r_energy_single);

    ##
    ## SIP
    ##

    Q_sip = mkmodel_sip( pp, aa, bb, [gamma_umts(i), gamma_wifi(j)], mu );
    q_sip = ctmc(Q_sip);
    tput_sip(i,j) = dot(q_sip, r_tput_sip);
    energy_sip(i,j) = dot(q_sip, r_energy_sip);

    ##
    ## ABPS
    ##

    Q_abps = mkmodel_abps( pp, aa, bb, [gamma_umts(i), gamma_wifi(j)], mu );
    q_abps = ctmc(Q_abps);
    tput_abps(i,j) = dot(q_abps, r_tput_abps);
    energy_abps(i,j) = dot(q_abps, r_energy_abps);

    printf(".");
  endfor
  printf("\n");
endfor

## Compute battery consumption
Q_single = mkmodel_single(pp(1), aa(1), bb(1), gg(1), mu(1));
L = ctmc_exps(Q_single,0:T-1,[1 0 0 0]);
batt_single = (L*r_energy_single')';

Q_sip = mkmodel_sip( pp, aa, bb, gg, mu );
L = ctmc_exps(Q_sip,0:T-1,[1 0 0 0 00 0 0]);
batt_sip = (L*r_energy_sip')';

Q_abps = mkmodel_abps( pp, aa, bb, gg, mu );
q0 = zeros(1,16); q0(1) = 1;
L = ctmc_exps(Q_abps,0:T-1,q0);
batt_abps = (L*r_energy_abps')';

##############################################################################
## THROUGHPUT
makeplot(gamma_wifi, gamma_umts, tput_single, tput_sip, tput_abps, "Throughput");
print("throughput.eps","-depsc2");

##############################################################################
## POWER CONSUMPTION
ax = makeplot(gamma_wifi, gamma_umts, energy_single, energy_sip, energy_abps, "Power Consumption (W)");
print("energy.eps","-depsc2","-tight");

makeplot_slices(gamma_wifi, gamma_umts, energy_single, energy_sip, energy_abps, ax,"energy_slices");

##############################################################################
## THROUGHPUT / POWER CONSUMPTION
ax = makeplot(gamma_wifi, gamma_umts, \
	      tput_single ./ energy_single, \
	      tput_sip ./ energy_sip, \
	      tput_abps ./ energy_abps, \
	      "Throughput / Power Consumption");
print("throughput_energy.eps","-depsc2");

makeplot_slices(gamma_wifi, gamma_umts, tput_single ./ energy_single, tput_sip ./ energy_sip, tput_abps ./ energy_abps, ax, "throughput_energy_slices");

##############################################################################
## BATTERY USAGE
figure(1);
clf;

## Global plot parameters
set(1,"paperorientation","landscape");
papersize=[5 3.5] * 0.6; 
margin=[0 0];
set(1,"papersize",papersize);
set(1,"paperposition", [margin papersize-margin*2]);

## interpolate
Elim = fix( min( [max(batt_single) max(batt_sip) max(batt_abps) ]) );

x_single = 1:Elim;
y_single = interp1(batt_single,0:T-1,x_single);

x_sip = 1:Elim;
y_sip = interp1(batt_sip,0:T-1,x_sip);

x_abps = 1:Elim;
y_abps = interp1(batt_abps,0:T-1,x_abps);


best_single = max( find( y_single >= y_abps & y_single >= y_sip) );
best_sip = max( find( y_sip >= y_single & y_sip >= y_abps ) );

plot(x_single,y_single,";UMTS;", "linewidth", 2, \
     x_sip, y_sip, ";SIP;", "linewidth", 2, \
     x_abps, y_abps, ";ABPS;", "linewidth", 2);
text(x_abps(100),y_abps(100),"ABPS \\rightarrow","horizontalalignment","right");
text(x_sip(50),y_sip(50),"SIP \\rightarrow","horizontalalignment","right");
ax = axis();
line( [best_single best_single], [0 ax(4)], "linestyle", "-" );
line( [best_sip best_sip], [0 ax(4)], "linestyle", "-" );
text(best_single, 5," x'","horizontalalignment","left");
text(best_sip, 5," x''","horizontalalignment","left");
axis("tic");
xlabel("Initial Capacity (x)");
ylabel("Battery duration ( R^{-1}(x) )");
legend("location","southeast");
print("battery_duration.eps","-deps2","-tight");
