## This script produces the plot for the case study described in the
## paper:
##
## M. Marzolla, R. Mirandola, Performance-Aware Reconfiguration of
## Software Systems, Technical Report UBLCS-2010-11, Department of
## Computer Science, University of Bologna, May 2010, available at
## http://www.cs.unibo.it/pub/TR/UBLCS/ABSTRACTS/2010.bib?ncstrl.cabernet//BOLOGNA-UBLCS-2010-11
##
## This script has been tested with GNU Octave version 3.2.4. You can
## download GNU Octave from http://www.octave.org/. This script also
## requires the qnetworks toolbox version 0.8.3, which can be downloaded
## from http://www.moreno.marzolla.name/software/qnetworks/
##
## This is version 1.0 of the script, written by Moreno Marzolla,
## <marzolla (at) cs.unibo.it> on march-may 2010.
##
## This script is released by the author in the public domain.

1; # this is not a function file

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

page_screen_output(0); # avoid output pagination

##############################################################################
## Given a square matrix SS and a vector idx, returns the vector S such
## that S(i) = SS(i,idx(i)). See
## https://www-old.cae.wisc.edu/pipermail/help-octave/2010-March/018479.html
function S = colidx(SS,idx)
  assert( rows(SS) == length(idx) );
  assert( all( idx >= 1 ) );
  assert( all( idx <= columns(SS) ) );
  S = SS(sub2ind( size(SS), 1:length(idx), idx ) );
endfunction

##############################################################################
## This function performs fractal interpolation of an array of integers
## values. Given a vector N, this function interpolate all values
## between N(i) and N(j) (j>i). Returns the new vector with values
## 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 graph. 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().
## Usage 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

##############################################################################
## Estimate the response time for model "model" with N requests, subject
## to configuration "conf". The response time is estimated using
## Balanced System bounds
function Rs = response_time_bounds( N, model, conf )
  [Xl Xu Rl Ru] = qnclosedbsb(N, colidx(model.S, conf));
  Rs = (Rl + Ru)/2;
endfunction

##############################################################################
## Compute the system response time using the QN model and MVA algorithm.
function Rs = response_time_mva( N, model, conf )
  [U R] = qnclosed(N, colidx(model.S, conf), ones(size(conf)));
  Rs = sum(R);
endfunction

##############################################################################
## Downgrade the component for which the ratio of service demand and
## utility is the maximum. This means that we prefer the component
## combining high service demand and low utility.
## 
## INPUTS
## N:           number of currently active requests
## model:       queueing network model as returned by the mkmodel function
## Rsmax:       maximum response time
function conf = downgrade( N, model, conf, Rsmax, resp_time_estimator )
  idx = find( conf>1 ); # idx = idx of components which can be degraded
  while( length(idx) > 0 )
    ## This is going to be tricky, first we need to project all data to
    ## the index set idx
    conf_proj = conf(idx); # configuration
    D_proj = model.S(idx,:); # service demand
    util_proj = model.util(idx,:);
    ## get maximum
    [U maxidx] = max(colidx(D_proj,conf_proj) ./ \
		     colidx(util_proj, conf_proj));
    B = idx(maxidx);
    ## downgrade component k
    conf(B)--;
    Rs = resp_time_estimator( N, model, conf );
    if ( Rs < Rsmax )
      break; # end reconfiguration
    else
      idx = find( conf>1 );
    endif
  endwhile
endfunction

##############################################################################
## Try to upgrade the component with minimum ratio D(k)/Util(k). D is the
## service demand of the upgraded component, Util is the utility of the
## upgraded component. Thus we prefer to upgrade components with the
## greatest utility and lowest demand.
##
## INPUTS
## N:           number of currently active requests
## model:       queueing network model as returned by the mkmodel function
## Rsmax:       maximum response time
function conf = upgrade( N, model, conf, Rsmax, resp_time_estimator )
  idx = find( conf<model.nvar ); # idx = idx of components which can be upgraded
  while( length(idx) > 0 )
    ## First we need to project all data to the index set idx
    conf_proj = conf(idx)+1; # new configuration
    D_proj = model.S(idx,:); # service demand
    util_proj = model.util(idx,:);
    ## minidx is the index such that D_proj(minidx)/util_proj(minidx) is minimum
    [U minidx] = min(colidx(D_proj,conf_proj) ./ \
		     colidx(util_proj, conf_proj));
    U = idx(minidx); # U is the component to upgrade
    assert( conf(U)<model.nvar(U) ); # sanity check
    ## ok, we have a candidate U. Try to upgrade it
    newconf = conf; newconf(U)++;
    Rs = resp_time_estimator( N, model, newconf );
    if ( Rs < Rsmax )
      conf = newconf;
      idx = find(conf<model.nvar);
    else
      ## upgrading U resulted in a violation of the response time
      ## limit. Thus, remove U from the candidate set, and try with another
      ## candidate.
      idx = idx(idx != U);
    endif
  endwhile
endfunction

##############################################################################
## Given a population size N and a model, check whether the maximum
## response time Rsmax is violated. If the system response time, as
## computed solving a QN model, is greater than Rsmax, then the
## downgrade() function is invoked. Otherwise, the upgrade() function is
## called.
##
## resp_time_estimator is a handle to a function which takes three
## parameters (N, model, conf) and returns a single scalar Rs which is
## the estimated response time for the model with N requests and given
## configuration. Currently, resp_time_estimator can be set either to
## response_time_mva or response_time_bounds.
##
## Returns the new configuration
function conf = adapt( N, model, conf, Rsmax, resp_time_estimator )
  ( all( conf <= model.nvar ) && all( conf >= 1 ) ) || \
      error("The current model configuration is not valid");
  (N > 0) || \
      error("The population size N must be >0");
  nargin == 5 || \
      print_usage();
  ## Compute the system response time for the current configuration
  Rs = response_time_mva(N, model, conf);
  if ( Rs > Rsmax )
    conf = downgrade(N, model, conf, Rsmax, resp_time_estimator );
  else
    conf = upgrade(N, model, conf, Rsmax, resp_time_estimator );
  endif
  ## disp(conf); # show the configuration
endfunction

##############################################################################
## Compute the utility of a specific configuration
function result = utility( util, conf )
  result = sum( colidx( util, conf ) );
endfunction

##############################################################################
## Build the model
function model = mkmodel( S, util, nvar )
  all( all( S(:,1:end-1) <= S(:,2:end) ) ) || \
    error("S() is not monotonically strictly decreasing" );
  all( all( util(:,1:end-1) <= util(:,2:end) ) ) || \
    error("util() is not monotonically strictly increasing" );
  model.S = S;
  model.util = util;
  model.nvar = nvar;
endfunction

##############################################################################
## Build a random model
##
## INPUTS:
## ncomp: number of software components
## nlev: number of quality levels offered by each component
##
## OUTPUT:
## model: software model as produced by mkmodel()
function model = mkrandommodel( ncomp, nlev )
  nvar = nlev*ones(1,ncomp);
  S = util = zeros(ncomp,nlev);
  S(:,1) = rand(ncomp,1)*10;
  util(:,1) = rand(ncomp,1)*10;
  for i=1:nlev-1
    S(:,i+1) = S(:,i) + rand(ncomp,1)*10;
    util(:,i+1) = util(:,i) + rand(ncomp,1)*10;
  endfor
  model = mkmodel( S, util, nvar );
endfunction

##############################################################################
## Model solution
##
## Given an input model, this function simulates adaptation with MVA and
## performance bounds, and computes the total utility in each case, and
## the sum of response time violations. The latter is defined as follows:
## let R(i) be the system response time at time i. The sum of response time
## violations is defined as sum( R( R>Rmax ) - Rmax ).
##
## INPUTS
## Rmax     = maximum allowed response time; if -1, then autocompute Rmax 
## model    = system model as created by mkmodel()
## N        = N(i) is the population size at time i
## plotname = name of the plot file; if omitted, no plot is produced
##
## Note that Rmax is autocomputed as the maximum response time which
## can be achieved, over the given population vector N, by setting all
## components at _worst_ quality level.
##
## OUTPUTS
## util_mva     = total utility with adaptation, MVA
## R_mva        = sum of response time violations, MVA
## util_pb      = total utility with adaptation, performance bounds
## R_pb         = sum of response time violations, performance bounds
## util_noadapt = total utility without adaptation, max quality
## R_noadapt    = sum of response time violations, no adaptation, max quality
## t_ave_mva    = avg time to reconfigure using mva
## t_ave_pb     = avg time to reconfigure using performance bounds
function [util_mva R_mva util_pb R_pb util_noadapt R_noadapt t_ave_mva t_ave_pb] = simulate( Rmax, model, N, plotname )
  idx = 1:length(N);
  if ( 2 == exist([plotname ".dat"]) )
    disp(["Reading data from file " plotname ".dat"] );
    load( [plotname ".dat"] );
  else
    Rs_pb = Rs_mva = Rs_noadapt = ut_pb = ut_mva = ut_noadapt = zeros(size(N));
    time_mva = time_pb = zeros(size(N));
    adapt_points = [];
    ## initial configuration is best quality level
    conf_mva = conf_pb = conf_best = model.nvar; 
    conf_worst = ones(size(model.nvar));
    
    R_best_min = response_time_mva(min(N), model, conf_best);
    R_best_max = response_time_mva(max(N), model, conf_best);
    R_worst_min = response_time_mva(min(N), model, conf_worst);
    R_worst_max = response_time_mva(max(N), model, conf_worst);
    
    printf("Response time @ best quality (min/max): %g/%g\n", 
           R_best_min, R_best_max );
    printf("Response time @ worst quality (min/max): %g/%g\n", 
           R_worst_min, R_worst_max );
    
    if ( Rmax < 0 ) 
      ## autocompute Rmax
      Rmax = R_worst_max;
    endif

    for i=idx
      ##
      ## Solve with MVA
      ##
      tic();
      conf_new_mva = adapt(N(i), model, conf_mva, Rmax, @response_time_mva);
      time_mva(i) = toc();
      if ( any(conf_new_mva!=conf_mva) )
        adapt_points = [adapt_points i];
        conf_mva = conf_new_mva;
      endif
      Rs_mva(i) = response_time_mva(N(i), model, conf_mva ); 
      ut_mva(i) = utility( model.util, conf_mva );
      
      ##
      ## Solve with performance bounds
      ##
      tic();
      conf_pb = adapt(N(i), model, conf_pb, Rmax, @response_time_bounds);
      time_pb(i) = toc();
      Rs_pb(i) = response_time_mva(N(i), model, conf_pb );
      ut_pb(i) = utility( model.util, conf_pb);
      
      ##
      ## Compute parameters without adaptation
      ##
      Rs_noadapt(i) = response_time_mva(N(i), model, conf_best);
      ut_noadapt(i) = utility( model.util, conf_best );
    endfor
    save([plotname ".dat"], "Rmax", "adapt_points", \
         "ut_mva", "Rs_mva", \
         "ut_pb", "Rs_pb", \
         "Rs_noadapt", "ut_noadapt", \
         "time_mva", "time_pb");
  endif

  util_mva = sum(ut_mva);
  R_mva = sum( Rs_mva( Rs_mva > Rmax ) - Rmax );
  util_pb = sum(ut_pb);
  R_pb = sum( Rs_pb( Rs_pb > Rmax ) - Rmax );
  util_noadapt = sum(ut_noadapt);
  R_noadapt = sum( Rs_noadapt( Rs_noadapt > Rmax ) - Rmax );
  t_ave_mva = mean(time_mva);
  t_ave_pb = mean(time_pb);

  printf("Total utility (MVA): %d\n", util_mva);
  printf("Total utility (bounds): %d\n", util_pb);
  printf("Integrated response time overflow (MVA): %f\n", R_mva );
  printf("Integrated response time overflow (bounds): %f\n", R_pb );

  if ( nargin == 4 ) 
    ## make plot
    set(gcf,"paperorientation","landscape");
    papersize=[5 4] * 0.9; margin=[0 0];
    set(gcf,"papersize",papersize);
    set(gcf,"paperposition", [margin papersize-margin*2]);
    
    ## top part: response time plot
    subplot(4,1,[1:2]);
    plot(idx, Rs_noadapt,";Static conf.;", "linewidth", 2, \
         idx, Rs_mva,";PARSY+MVA;", "linewidth", 2, \
         idx, Rs_pb, ";PARSY+BSB;", "linewidth", 2, \
         adapt_points, Rs_mva(adapt_points), "+;Reconf. point;", "markersize", 3, "linewidth", 2, \
         idx, Rmax*ones(size(idx)), ";R_{max};");
    axis("labely");
    ylabel("System Response Time");
    ax = axis();
    ax(3) = 0;
    axis(ax);
    
    ## center part: utility
    subplot(4,1,3);
    plot(idx,ut_noadapt, ";Static conf.;", "linewidth", 2, \
         idx,ut_mva, ";PARSY+MVA;", "linewidth", 2, \
         idx,ut_pb, ";PARSY+BSB;", "linewidth", 2 \
         );
    ylabel("Utility");
    legend("location","southeast");
    ax = axis();
    ax(3) = 0;
    axis(ax);
    axis("labely");
    
    ## bottom part: population size
    subplot(4,1,4);
    plot(idx,N, "linewidth", 2);
    xlabel("Time step");
    ylabel("N. Active Users");
    
    print([plotname ".eps"], "-deps2", "-mono");
  endif
endfunction

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

##
## Create table
##
table = fopen("ex_table.tex", "w");

N = fractal_build_population([  1 20 40 60 80 100 120 140 160 180 200],
			     [ 30 20 10 80 100 50  20 40  60  10 20] );

i = 1;
ut = rr = zeros([12 3]);
lab = cell(1,12);
Rmax = -1;
for ncomp=[5 10 20 50]
  for nlev=[2 3 5]
    printf("***\n*** %d components, %d levels each\n***\n", ncomp, nlev);
    model = mkrandommodel( ncomp, nlev );
    [u_mva R_mva u_pb R_pb u_static R_static time_ave_mva time_ave_pb] = simulate( Rmax, model, N, sprintf("example-%d-%d", ncomp, nlev) );
    ut(i,:) = [u_mva u_pb u_static];
    rr(i,:) = [R_mva R_pb R_static];
    lab{i} = sprintf("%d/%d",ncomp,nlev);
    i++;
    fprintf(table, "%d & %d & %.2f & %.2f & %.2f & %.2f & %.2f & %.2f\\\\\n", 
            ncomp, nlev, u_static, R_static, u_mva, R_mva, u_pb, R_pb);
  endfor
endfor

set(gcf,"paperorientation","landscape");
papersize=[7 4] * 0.75; margin=[0 0];
set(gcf,"papersize",papersize);
set(gcf,"paperposition", [margin papersize-margin*2]);
subplot(2,1,1)
bar(ut);
title("Utility");
h = gca();
set(h,"xtick",1:12);
set(h,"xticklabel",lab);
legend({"PARSY+MVA","PARSY+BSB","Static"}, "location", "northwest");

subplot(2,1,2);
bar(rr);
title("Response time overflow (log scale)");
h = gca();
set(h,"yscale","log");
set(h, "xtick",1:12);
set(h, "xticklabel",lab);
legend({"PARSY+MVA","PARSY+BSB","Static"}, "location", "northwest");
print("hist.eps","-color");