## This script implements the dynamic reconfiguration algorithm
## described in the paper: "Dynamic Scalability for Next Generation
## Gaming Infrastructures", by G. d'Angelo, S. Ferretti, M. Marzolla,
## Technical Report UBLCS-2010-12, Dipartimento di Scienze dell'Informazione,
## Universita' di Bologna, Italy.
## Specifically, this script has been used to produce the results
## described in the paper above.
##
## This script has been tested with GNU Octave version 3.2.4.
## You can download GNU Octave from http://www.octave.org/. 
##
## This script was written by Moreno Marzolla,
## <moreno.marzolla (at) unibo.it> on october--december 2010.
##
## This script is released by the author according to the Creative
## Commons CC0 1.0 Universal license
## http://creativecommons.org/publicdomain/zero/1.0/

##############################################################################
##
## NOTE: As explained in the paper, in order to reduce the number of
## free parameters, in this simulation we assume that the number of
## concurrent users in the system is the only parameter. We set the
## average service time at each service center, and assume that each
## request visits each tier exactly once. In this way we can compute
## everything needed by the reconfiguration algorithm as a function of
## the number of concurrent users only.
##
##############################################################################

1; # this is not a function file

figure (1, "visible", "off"); # do not show plot window

page_screen_output(0); # avoid output pagination

##############################################################################
## Mean Value Analysis (MVA) for closed, product-form single-class
## queueing networks. The computational cost of this function is O(
## N*length(D) ).
##
## Inputs:
## N    Number of requests in the network (scalar)
## D    D(k) is the service demand at center k
##
## Outputs:
## U    U(k) is the estimated utilization of center k
## Rs   (scalar) estimated residence time at center k
function [U Rs] = mva(N, D)
  Q = zeros(size(D));
  for n=1:N
    R = D .* (1 + Q);
    X = n / sum(R);
    Q = X*R;
  endfor
  U = X*D;
  Rs = sum(R);
endfunction

##############################################################################
## Balanced System Bound (BSB) estimates for closed, product-form
## single-class queueing networks. The computational cost of this
## function is O( length(D) )
##
## Inputs:
## N    Number of requests in the network (scalar)
## D    D(k) is the service demand at center k
##
## Outputs:
## U    U(k) is the _approximate_ utilization of center k
## Rs   (scalar) _approximate_ residence time at center k
function [U Rs] = bsb(N, D)
  D_max = max(D);
  D_tot = sum(D);
  D_ave = mean(D);
  Xl = N/(D_tot+( (N-1)*D_max ) );
  Xu = min( 1/D_max, N/( D_tot+( (N-1)*D_ave ) ) );
  Rl = max( N*D_max, D_tot+( (N-1)*D_ave ) );
  Ru = D_tot + (N-1)*D_max;
  X = mean([Xl Xu]);
  Rs = mean([Rl Ru]); 
  U = X*D;
endfunction

##############################################################################
## Compute the response time of the cloud-based gaming system using
## a QN model. 
##
## Inputs:
## N		(scalar) number of concurrent users
## conf         (vector) [num_gateways num_cell_srv num_db_srv]
## model        (struct) structure as returned by mkmodel
## mode         (string) if "exact", evaluate the model using MVA (slower); if "approx" evaluate the model using BSB (faster, less precise). Default "approx"
##
## Outputs:
## U            U(k) is the utilization of tier k (k=1:3)
## Rs           (scalar) system response time
function [U Rs] = analyze_system(N, conf, model, mode );
  if ( nargin == 3 ) 
    mode = "approx";
  endif
  S = [model.SG model.SC model.SDB];
  V = 1 ./ conf;
  if ( strcmp(mode,"exact") ) 
    [U Rs] = mva(N, S.*V);
  elseif ( strcmp(mode, "approx") )
    [U Rs] = bsb(N, S.*V);
  else
    error("Invalid mode %s: must be \"exact\" or \"approx\"", mode);
  endif
endfunction

##############################################################################
## Acquire resources
function newconf = acquire( N, conf, model );
  newconf = conf;
  thr = mean([model.Rlow model.Rhigh]); # threshold
  do
    [U Rsys] = analyze_system(N, newconf, model );
    if ( Rsys >= thr )
      [Umax Umaxidx] = max(U); # identify bottleneck
      newconf(Umaxidx)++;
    endif
  until (Rsys < thr);
endfunction

##############################################################################
## Release resources
function newconf = release( N, conf, model )
  candidates = find(conf > 1);
  thr = mean([model.Rlow model.Rhigh]); # threshold

  while( length(candidates) > 0 )
    G = conf(1);
    C = conf(2);
    DB = conf(3);

    [U Rsys] = analyze_system(N, conf, model );

    [Umin minidx] = min( U(candidates) );
    Uminidx = candidates(minidx);

    newconf = conf; newconf(Uminidx)--;
    [U Rsys] = analyze_system(N, newconf, model );
    if ( Rsys < thr )
      conf = newconf;
      candidates = find(conf > 1);
    else
      candidates = candidates( candidates != Uminidx );
    endif
  endwhile
  newconf = conf;
endfunction

##############################################################################
## Given a vector N of integers, this function interpolate all values
## between N(i) and N(j) (j>i). Returns the new vector with the same length
## of N, with all values in the range [i,j] interpolated.
##
## Example:
##
## N = zeros(1,10);
## N(1) = 10; N(10) = 5;
## N = fractal_interpolate(N,1,10);
function N = fractal_interpolate(N, i, j)
  if ( j > i+1 )
    m = round((i+j)/2);
    d = max(abs(N(i)-N(j)),2);
    v = round((N(i)+N(j))/2 + unifrnd(-d/2, d/2));
    N(m) = v;
    N = fractal_interpolate(N,i,m);
    N = fractal_interpolate(N,m,j);
  endif
endfunction

##############################################################################
## Produce a fractal population sequence. N(idx) = val is the initial value
## (both idx and val must be vectors of the same size). All missing
## values are interpolated using the function fractal_interpolate().
##
## Example
## 
## fractal_build_population([1 10 20], [10 30 5])
##
## This sets N(1) = 10, N(10) = 30 and N(20) = 5; all remaining values
## are interpolated.
function N = fractal_build_population(idx, val)
  N = 0;
  N(idx) = val; # initial values
  for i=2:length(idx)
    N = fractal_interpolate(N,idx(i-1),idx(i));
  endfor
endfunction

##############################################################################
## Compute the moving average on vector X as
## ( X(i) + X(i-1) + X(i-2) + ... + X(i-W+1) ) / W
## if i<W, the compute
## ( X(i) + X(i-1) + ... + X(1) ) / i
function result = movavg( X, i, W)
  result = mean( X( max(1,i-W+1) : i ) );
endfunction

##############################################################################
## Build a structure containing all parameters for the datacenter model
##
## Inputs:
## Rmax:  Max response time
## Rhigh: High threshold 
## Rlow:  Low threshold
## SG:    service times at gateways
## SC:    service times at cell servers
## SDB:   service times at DB servers
##
## Outputs:
## model  Data structure containing the above information
function model = mkmodel( Rmax, Rhigh, Rlow, SG, SC, SDB )
  assert( Rmax > SG + SC + SDB );
  model.Rmax = Rmax;
  model.Rhigh = Rhigh;
  model.Rlow = Rlow;
  model.SG = SG;
  model.SC = SC;
  model.SDB = SDB;
endfunction

##############################################################################
## Perform a simulation run. This function produces three filds in the local
## directory:
## - <plotname>.dat: this file contains the raw data produced by the simulation 
##   run. NOTE: IF THIS FILE ALREADY EXISTS WHEN THIS FUNCTION IS CALLED, THEN
##   ALL SIMULATION RESULTS ARE READ FROM THE FILE INSTEAD OF BEING RECOMPUTED.
##   MAKE SURE TO REMOVE <plotname>.dat IF YOU WANT TO RE-RUN THE SIMULATION
##   WITH DIFFERENT PARAMETERS.
## - <plotname>_color.eps: graph showing the results.
##
## Inputs:
## N		N(i) is the number of users at step i
## conf		conf(i) is the number of nodes initially allocated to tier i
## W		window length used for computing moving averages
## SG		Average service time at gateways
## SC		Average service time at cell servers
## SDB		Average service time at DB servers
## Rmax		Max allowed response time
## Rhigh	"hard" threshold
## Rlow		"soft" threshold
## plotname	basename of the figures to produce
##
## Outputs:
## nv           Number of violations of the constraint R<R_max
## ro           Response time overflow, defined as the Area between the 
##              line y=R_max and the response time y=RS(t)
## nr           Number of reconfigurations
## ts           Total number of servers allocated during the whole period
## lmean        Average duration of response time violation intervals
## lmax         Max duration of response time violation intervals
function [nv ro nr ts lmean lmax] = simulate( N, conf, W, SG, SC, SDB, Rmax, Rhigh, Rlow, plotname )
  idx = 1:length(N);
  if ( 2 == exist([plotname ".dat"]) )
    disp(["Reading data from file " plotname ".dat"] );
    load( [plotname ".dat"] );
  else
    Rs = Ravg = zeros(size(N));
    adapt_points = [];
    confs = zeros(length(N),3); # confs(i,:) is the configuration at time i

    ## create a model structure
    model = mkmodel( Rmax, Rhigh, Rlow, SG, SC, SDB );

    printf("Legend (fields listed from left to right)\n\n");
    printf("(1) Index (timestamp)\n");
    printf("(2) Number of concurrent requests\n");
    printf("(3) Current Response Time\n");
    printf("(4) Time-averaged Response Time\n");
    printf("(5) Acquire / Release resources\n");
    printf("(6) New configuration (G, C, DB)\n\n");
    for i=idx
      [U Rs(i)] = analyze_system( N(i), conf, model, "exact");
      confs(i,:) = conf(:);

      ## Perform reconfiguration if necessary
      Ravg(i) = movavg(Rs,i,W);
      action = "none";
      if ( rem(i,W) == 1 ) 
        if ( Ravg(i) < Rlow )
	  conf = release(N(i), conf, model );
          adapt_points = [adapt_points i];
          action = "release";
        elseif ( Ravg(i) > Rhigh )
	  conf = acquire(N(i), conf, model );
          adapt_points = [adapt_points i];
          action = "acquire";
        endif
      endif

      printf("%3d\t%6d\t%9.3f\t%9.3f\t%7s\t(%2d, %2d, %2d)\n",
             i, N(i), Rs(i), Ravg(i), action, conf(1), conf(2), conf(3) );
    endfor
    ## Dump all simulation results to the file <plotname>.dat
    save([plotname ".dat"], \
	 "Rmax", "N", "adapt_points", "Rs", "confs", "Ravg", "Rmax", "Rhigh", "Rlow");
  endif

  nv = sum(Rs > Rmax);
  ro = sum(Rs( Rs > Rmax ) - Rmax);
  nr = length(adapt_points);
  ts = sum(sum(conf));

  ll = [];
  Roverflow = Rs>Rmax; # boolean vector
  cur_len = 0;
  for i=1:length(Rs)
    if (Roverflow(i) == 1)
      cur_len++;
    else
      if ( cur_len > 0 )
        ll = [ll cur_len];
      endif
      cur_len = 0;
    endif
  endfor
  if ( cur_len > 0 )
    ll = [ll cur_len];
  endif

  lmean = mean(ll);
  lmax = max(ll);
  printf("\nSummary:\n\n");
  printf("Response time overflow    : %f\n", ro );
  printf("Avg. violation duration   : %f/%f\n", lmean, lmax );
  printf("Max violation %%           : %f\n", (max(Rs(Rs > Rmax)) - Rmax)/Rmax );
  printf("Number of violations      : %d/%d\n", nv, length(N) );
  printf("Number of reconfigurations: %d\n", nr);
  printf("Total allocated servers   : %d\n", ts);

  ##
  ## Produce the plot
  ##
  set(gcf,"paperorientation","landscape");
  papersize=[5 4] * 0.85; 
  ##papersize=[5 4] * 1.2;
  margin=[0 0];
  set(gcf,"papersize",papersize);
  set(gcf,"paperposition", [margin papersize-margin*2]);
  
  ## top part: response time
  subplot(10,1,[1:5]);
  plot(idx, Rs,";Response Time;", "linewidth", 2, \
       idx, Ravg, ";Moving Av. Resp. Time;", "linewidth", 1.5, \
       adapt_points, Rs(adapt_points), "+;Reconf. point;", "markersize", 3, "linewidth", 2,\
       idx, Rmax*ones(size(idx)), "k", "linewidth", 1.5, \
       idx, Rhigh*ones(size(idx)), "k", \
       idx, Rlow*ones(size(idx)), "k"
       );
  legend("location", "southeast");
  axis("labely");
  ylabel("Response Time");
  text(max(idx),Rmax," R_{max}");
  text(max(idx),Rhigh, " R_{high}");
  text(max(idx),Rlow, " R_{low}");
  ax = axis();
  ax(3) = 0;
  ax(4) = Rmax*1.4;
  axis(ax);
  
  ## middle part: number of instances of each type
  subplot(10,1,[6 7]);
  h = area(idx, confs);
  set(h(1),"FaceColor",0.4*[1 1 1]);
  set(h(2),"FaceColor",0.6*[1 1 1]);
  set(h(3),"FaceColor",0.8*[1 1 1]);
  axis("labely");
  ylabel("C_1, C_2, C_3");
  ax = axis();
  ax(3) = 0;
  ##ax(4) = 1.5*max(sum(confs,2));
  axis(ax);
  ##legend({"C_1", "C_2", "C_3"},
  ##       "location", "northwest"); #, "orientation", "horizontal");

  ## bottom part: population size and delta
  subplot(10,1,[8 9]);
  plot(idx, N, "linewidth", 2);
  axis("labely");
  ylabel("N. Active Users");
  ax = axis();
  ax(3) = 0;
  axis(ax);

  ## plot delta; each bar in this graph is the difference between
  ## adjacent population sizes
  subplot(10,1,10);
  stem(idx,[0 diff(N)],"markersize",0,"linewidth",3);
  axis("ticx");
  xlabel("Time step");

  ## print([plotname "_mono.eps"], "-deps2", "-mono" );
  print([plotname "_color.eps"], "-deps2", "-color");
endfunction

##############################################################################
## initialize random number generator
rand("state",[13 9 -11 71 681653 2 9 4 1 9 10 7 -23 8762 467256 481 1 0] );

fout = fopen("results.tex", "w");

##############################################################################
##############################################################################
###
### Example 1
###
##############################################################################
##############################################################################

printf("***\n*** EXAMPLE 1\n***\n");

## Build the population 
N = fractal_build_population([   1   15   35   50   80  90  120 140  150  180  200],
			     [ 300  500  650  400  650  800 650 300  400  650  400]);

[nv ro nr ts lmean lmax] = simulate( N, [5, 5, 5], 5, 0.08, 0.8, 0.58, 100, 90, 70, "example1" );
fprintf(fout,"Test~\\ref{fig:ex1} & %d--%d & %d & %d & %.2f/%.2f\\\\\n", min(N), max(N), nr, nv, lmean, lmax);

x = [1:10:200 200];

##############################################################################
##############################################################################
###
### Example 2
###
##############################################################################
##############################################################################

printf("***\n*** EXAMPLE 2\n***\n");
y = 10000 + sin(x/20) .* 7000 + 1990*rand(1,length(x));
N = fractal_build_population(x, y);

[nv ro nr ts lmean lmax] = simulate( N, [10, 120, 90], 5, 0.08, 0.8, 0.58, 100, 90, 70, "example2" );
fprintf(fout,"Test~\\ref{fig:ex2} & %d--%d & %d & %d & %.2f/%.2f\\\\\n", min(N), max(N), nr, nv, lmean, lmax);

##############################################################################
##############################################################################
###
### Example 3
###
##############################################################################
##############################################################################

printf("***\n*** EXAMPLE 3\n***\n");
y = 10000 + sin(x/15) .* 7000 + 1990*rand(1,length(x));
N = fractal_build_population(x, y);

[nv ro nr ts lmean lmax] = simulate( N, [10, 120, 90], 5, 0.08, 0.8, 0.58, 100, 90, 70, "example3" );
fprintf(fout,"Test~\\ref{fig:ex3} & %d--%d & %d & %d & %.2f/%.2f\\\\\n", min(N), max(N), nr, nv, lmean, lmax);

##############################################################################
##############################################################################
###
### Example 4
###
##############################################################################
##############################################################################

printf("***\n*** EXAMPLE 4\n***\n");
y = 10000 + sin(x/10) .* 7000 + 1990*rand(1,length(x));
N = fractal_build_population(x, y);

[nv ro nr ts lmean lmax] = simulate( N, [10, 140, 100], 5, 0.08, 0.8, 0.58, 100, 90, 70, "example4" );
fprintf(fout,"Test~\\ref{fig:ex4} & %d--%d & %d & %d & %.2f/%.2f\\\\\n", min(N), max(N), nr, nv, lmean, lmax);

fclose(fout);

