%Data inserted for a layered queue model of a web application, works with
%modelCalc_LQWeb
% Feb 5 2010 regressModel adapted from version for WOSP/SIPEW tutorial
% Generalized and formalized, and input parameter selector introduced.

%Jan 18 2010 regressModel adapted from regressrun_lqn_exp3_all_dual.m
%calls modelSens for model results for 2 Queue model
%template [h,H]=modelSens(param,K,C,Npop,selectpar)
%h = values, H = sensitivities over selected params

% enter inputs as follows:
% Model-related:
% (1) nPotential = integer number of parameters passed to modelCalc, including fitted, controlled and fixed.
%     = size(param) = size(parInit) = size(parNames) = size(parFit) = size(parControlled) = size(parInteger)
% (2) parNamesCell = ROW cell vector of string names for the parameters
% (3) parInteger = col vector with 1 for integer parameters, else 0
% (4) parInit = initial parameter col vector, with values for fixed parameters
% (5) parFit = col vector designating fitted parameters by 1, otherwise 0
% (6) pControlInd = row vector connecting controlled parameters to data columns
%         (elements are the columns of the data, which are controlled
%         parameters... number of elements = number of controlled
%         parameters)
% 
% Data-related:
% (7) TMAX = the number of trials
% (8) mpotential = the number of data values per trial, including controlled parameters and measured data
%                 = size(measUsed)
% (9) dataNamesCell = a ROW cell vector of string names for the data values
% (10) dataUsed = a column vector of size mpotential designating just the data values 
% which are measured performance values used for fitting
%           dataUsed(i) = 1 iff data column i is used
% %The model calculation should return a vector h of size = mpotential, with
% %zeros for the controlled parameter positions and performance measures
% %corresponding by position to the measured values.
% (11) data = a TMAX by mpotential array of data values. Each row is the data for one trial
%       CONTROLLED parameters should be the first elements of each row.
% (12) dataVariances = a mpotential-column-vector of relative variances of the data
% elements (with zeros for the controlled parameters)

% Program control parameters:
% (13) maxIter = integer limit on the regression iterations (suggested value about 20, make it bigger if it meets this limit)
% (14) convtest = a small convergence test value such as .00001
% Special data for this model
% In the example, K = no of queues, C = classes

% Now put the data. Data for the two class example is given as a guide
%Model-related:
nPotential = 9;
parNamesCell = {'  $nusers', '  $DirDemand', '  $WSDemand', '  $appDph1', '  $appDph2', '  DBopD', '  mApp', '  $wsThreads','  $think'};  %parameter names in template_LQWeb.xlqn
parInteger = [1,0,0,0,0,0,1,1,0]';  %users, app and ws threads are integers
parInit = [10,1,2,1,1,1,2,1,5]';  %10 and 5 are fixed (10 is also controlled)
%parexact = [10,0.1,0.4,0.2,0.2,0.1,1,3,5]';  %used in the simulation
parFit = [0,1,1,1,1,1,1,1,0]';  %all but users and think times are fitted
pControlInd = [1];  %in the data matrix, col 1 $nusers is controlled parameter = 5*trial number

% Data-related:
TMAX = 10; %trials (the number of rows of data
mpotential = 6; %data elements per trial, the number of columns of data
%dataNames = [ '$nusers', '*********$N2', '$T11', '$T21', '$T31', '$n11', '$n21', '$n31', '$U11', '$U21', '$U31', '$Thru1', '$T12', '$T22', '$T32', '$n12', '$n22', '$n32', '$U12', '$U22', '$U32','$Thru2']';
dataUsed = [0,1,1,1,1,1]';
% from simulation
data=[5,0.8339 0.887684 0.351432 0.735055 0.131008
10,1.515 1.44121 0.356516 1.77054 0.247736
15,1.9875 2.47965 0.359583 2.56935 0.312413
20,2.1413 4.2804 0.356744 2.92819 0.331579
25,2.1936 6.37935 0.355907 2.99356 0.332453
30,2.1679 9.01488 0.353192 2.99926 0.335884
35,2.1405 11.3755 0.362741 2.99997 0.340277
40,2.1455 13.5649 0.35758 2.99995 0.337287
45,2.1344 15.9101 0.36153 2.99973 0.339217
50,2.1747 17.9796 0.358798 2.99987 0.33733];

%the following variances were taken from the confidence intervals reported
%by the simulation (square of half conf int/1.96)
dataVariances = [0,.257,.718,.01,.5912,.01334]; %1st col is controlled, has no variance

% Iteration Control parameters
maxIter = 50;
convtest = 0.00001;

% Special for this model:
% global K
% global C
% K = 3; %no of queues
% C = 2; %no of classes

% Then checks on sizes
[s1,s2] = size(data);
if (s1 ~= TMAX) 
    fprintf(1,'wrong number of trials (data rows)   \n', s1);
end
if(s2 ~= mpotential) 
    fprintf(1,'wrong number of datavalues in a trial   \n', s2);
end
if(size(dataUsed) ~= mpotential)
    fprintf(1,'wrong number of data-used indicators   \n', size(dataUsed));
end
if (size(parInit)~=nPotential)
    fprintf(1,'wrong number of inital parameter values   \n', size(parInit));
end
if (size(parInteger)~=nPotential)
    fprintf(1,'wrong number of integer parameter indicators   \n', size(parInteger));
end
if (size(parFit) ~= nPotential)
    fprintf(1,'wrong number of fitted parameter indicators   \n', size(parFit));
end


% Then derivation of Spar, S
% Spar = selector matrix for parameters, Spar*param = vector of fitted
% parameters
n = sum(parFit);
parFitt = parFit;
Spar = zeros(n,nPotential);
i = 1;
for j = 1:nPotential
    if (parFitt(j)>0)
        Spar(i,j) = 1;
        i = i + 1;
    end
end

%S = selector matrix for data used, S*data(t,:) = vector of data used for
%trial t
dataUsedd = dataUsed;
m = sum(dataUsed);
S = zeros(m,mpotential);
i = 1;
for j = 1:mpotential
    if (dataUsed(j)>0)
        S(i,j) = 1;
        i = i + 1;
    end
end

% Set up the weighting matrix, inverse of the measurement variance
% matrix when restricted to selected measurement values
Rbase = diag(dataVariances);
sqrtRbase = sqrt(Rbase);
R = S*Rbase*S';    %restrict to selected measurements
W = eye(m)/R;    %weight matrix = 1/var
sqrtR = S*sqrtRbase*S';

%iteration initialize
M=m;
N = n;
Tstudent = 1.96;  %critical value for conf int at 95% level, large DF.
partrial=parInit;
x = Spar*parInit;   %initial fitting vector
parIntSelected = Spar*parInteger;
xx=zeros(2*TMAX,N);
xx(1,:)=x;  %copies the column vector into a row!
diff=1;
sseold=inf;
sse=0;
P = zeros(N,N);
Pold=zeros(N,N);
xold=zeros(N,1);
bpmold=zeros(N,1);
mseold=0;

%begin iteration
iter=0;
while ((diff>convtest)&(iter<maxIter))
    iter=iter+1
    if(iter>2)
        sseold=sse;
        Htempold=Htemp;
        mseold=mse;
        Pold=P;
        xold=x;
        bpmold=bpm;
    end;
    
    P = zeros(N,N);
    V = zeros(1,N)';
    sse = 0;
    h=zeros(TMAX,M);
    for t = 1:TMAX    %sum over measurement steps
        %set populations from data, first C colunms        
        
        %set up parameters that are estimated, in partrial
        jpar=0;
        for ipar=1:nPotential
            partrial(ipar) = parInit(ipar);
            if(parFit(ipar)>0)
                jpar=jpar+1;
                partrial(ipar)= x(jpar); %set the values to use, from the estimated vector
            end;
        end;
        
        controlled = data(t,pControlInd);
        [hh,HH] = modelSens(partrial,controlled,parFit,parInteger);
        %hh=hh';
        hold=(S*hh)';  %row vector
        h(t,:) = hold;       
        Htemp = (S*HH*Spar'); %select sensitivities of measurements wrto estimated parameters
        etemp = (S*data(t,:)' - h(t,:)'); %column vector
        P = P + Htemp'*W*Htemp;
        V = V + Htemp'*W*etemp;
        sse = sse + etemp'*W*etemp;
    end;    %sums over t
%    sse,P,V
    
    b = (eye(N)/P)*V;   %the increment to the parameter vector
    mse = (sse - b'*V)/(M*TMAX - N);
    %limited change in direction of zero, prevent negative x
    % limited to a factor of 5 increase or decrease.
    trial = x + b;
    minx = 0.2*x;
    maxx=5*x;
    for i=1:N
        if(trial(i)>minx(i))
            x(i)=trial(i);
        else
            x(i)=minx(i);
        end;
        if(x(i)>maxx(i))
            x(i)=maxx(i);
        end;
    end;
% round off the integer parameters
    for i=1:N
        if(parIntSelected(i)>0)
            x(i) = round(x(i));
            x(i) = max([1,x(i)]);
        end;
    end;
 
    xx(iter,:)=x;
    %P = covar matrix estimate for b
    P = mse*eye(N)/P;
    bpm = Tstudent*diag(sqrt(P));%corrected from (1/mse)
    %brange(iter,:) = bpm; % the half confidence interval
    diff = ones(1,N)*abs(b);
    jpar=0;
    for ipar=1:nPotential
       par(ipar) = parInit(ipar);
       if(parFit(ipar)>0)
               jpar=jpar+1;
               par(ipar)= x(jpar);
        end;
    end; 
%    par
end;   %iteration
%if stopped on worse result, restore previous values
if(sseold<sse)
    'stopped on step that got worse'
    sse=sseold;
    Htemp=Htempold;
    mse=mseold;
    P=Pold;
    x=xold;
    bpm=bpmold;
end;

x
count=fprintf(1,'converged at %d iterations, diff = %f \n', iter, diff );
count=fprintf(1, 'results for estimated parameters and 95 percent conf int plusminus \n');

parFitPosition = zeros(1,nPotential);
j = 1;
for i = 1:nPotential;
    if (parFit(i) == 1)
        parFitPosition(j) = i;
        j = j+1;
    end
end

%print results
%parFitPosition(1:(j-1))
%parNamesCell
parNamesFitCell = parNamesCell(parFitPosition(1:(j-1)));
%parNamesFitCell
parNamesFit = char(parNamesFitCell);  %make a character array for printing
parNamesFit = reshape(parNamesFit',1,prod(size(parNamesFit))); %makes it a horizontal string
count=fprintf(1,'\n paramValue');
count=fprintf(1, ' %s',parNamesFit, '  '); %the final space seems to be necessary to get the text on a single line
count=fprintf(1,'\n paramValue');
count=fprintf(1,' % 7.3f ',x);
count=fprintf(1,'\n confInt+- ');
count=fprintf(1,' % 7.3f ',bpm);
count=fprintf(1,'\n');

%DONE


%Example of contour plot of sse for two parameters D21, D22: variable suffixed Con

%Contour plot
%paramCon = [2,3,4,10,5,10,1,30]';  %exact values
% set up the grid
% for iCon=1:40
%     paramCon(2) = 2.5+0.025*iCon
%     for jCon = 1:40
%         paramCon(6) = 9+.05*jCon;
%         sse=0;
%         
%         for t = 1:TMAX    %sum over measurement steps
%             %set populations from data, first C colunms        
%     %         
%     %         %set up parameters that are estimated, in partrial
%     %         jpar=0;
%     %         for ipar=1:nPotential
%     %             partrial(ipar) = paramCon(ipar);
%     %             if(selectpar(ipar)>0)
%     %                 jpar=jpar+1;
%     %                 partrial(ipar)= x(jpar); %set the values to use, from the estimated vector
%     %             end;
%     %         end;
%     %         
%     %        [hh,HH] = lqnSolve(partrial,selectpar,nPotential,m);
%             [hh,HH] = modelSens(paramCon,K,C,Npopp(t,:),parFit);
%             %hh=hh';
%             hold=(S*hh)';  %row vector
%             h(t,:) = hold;       
%             %size(h(t,:))
%             %size(z(t,:))
%             %Htemp = (S*HH*Spar'); %select sensitivities of measurements wrto estimated parameters
%             etemp = (S*z(t,:)' - h(t,:)'); %column vector
%             %P = P + Htemp'*W*Htemp;
%             %V = V + Htemp'*W*etemp;
%             sse = sse + etemp'*W*etemp;
%             sseCon(iCon,jCon)=sse;
%         end;
%     end;
% end;
% contour(sseCon,20); %20 contour lines, matlab chooses the vales


% % Example of printing and Plotting the prediction confidence intervals and the data points
% %conf interval for predictions at given value of t (= N, the user population)
% % plot over trials to 2*TMAX... in this experiment that means
% % extrapolating to higher populations, beyond the data
% 'prediction conf intervals at different populations t'
% 'pop    par1  plusminus    par2   plusminus   par3   plusminus'
% respt = zeros(2*TMAX,6);
% for t = 1:2*TMAX
%     [hh,HH] = mvaproc_think (x,t);
%     hh = S*hh;
%     HH=S*HH;
%     C = HH*P*HH';
%     for i = 1:M
%         ypm(i) = Tstudent*sqrt(C(i,i));
%     end;
%     tmin = (min([t,TMAX]));  %use this for the data points, which only go to TMAX
% % respt(t,6) is the computed value of response time from the throughput
%     respt(t,:)=[t, t/hh(4)-10, t/(hh(4)+ypm(4))-10, t/(hh(4)-ypm(4))-10, tmin, tmin/data(tmin,10)-10];
%     %[ t, hh(1),ypm(1), hh(2),ypm(2), hh(3), ypm(3) hh(4) ypm(4)];
%     count=fprintf(1,'% 7.3f % 7.3f % 7.3f % 7.3f % 7.3f % 7.3f % 7.3f % 7.3f % 7.3f \n',t, hh(1),ypm(1), hh(2),ypm(2), hh(3), ypm(3), hh(4), ypm(4));
% end;
% 
% plotmax= 2*TMAX; 
% plot(respt(1:plotmax,1),respt(1:plotmax,2),'k+-',respt(1:plotmax,1),respt(1:plotmax,3),'k+--',respt(1:plotmax,1),respt(1:plotmax,4),'k+--',respt(1:plotmax,5),respt(1:plotmax,6),'rx');
% xlabel('Population');
% ylabel('Predicted Response Time');
% %title('Predicted Response Time with 95% Confidence Limts');
% % %solid line for the step value function, dotted for smooth
% print -djpeg  QRespPred_think-data3.jpg;
% 
