% 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:
%1Q Data input for an example in regressModel:
nPotential = 4;
parNamesCell = {'  D1  ', 'D2  ','  D3  ','  Z  '};
parInteger = [0,0,0,0]';
parInit = [1,1,1,10]';
parFit = [1,1,1,0]';
pControlInd = [1];
TMAX = 10;
mpotential = 11;
dataNamesCell = {'  N  ', '  T1  ', '  T2  ','  T3  ',' n1  ',' n2 ',' n3 ',' U1 ',' U2 ',' U3 ',' f '};             

data = [1 1.21384 4.73153 2.64049 0.060692 0.2365765 0.1320245 0.0594781 0.264966 0.16847 0.05
2 1.60529 2.32118 4.72039 0.17337132 0.25068744 0.50980212 0.14541 0.204289 0.473156 0.108
3 1.71797 4.4539 4.26163 0.25425956 0.6591772 0.63072124 0.210235 0.432874 0.538983 0.148
4 2.82288 5.48302 6.825 0.4093176 0.7950379 0.989625 0.303433 0.586554 0.658857 0.145
5 3.49714 6.47212 8.42731 0.60850236 1.12614888 1.46635194 0.372242 0.499684 0.788018 0.174
6 3.00036 4.3094 13.6851 0.6150738 0.883427 2.8054455 0.389882 0.523607 0.875007 0.205
7 3.0938 5.55339 11.3586 0.7486996 1.34392038 2.7487812 0.439166 0.599545 0.870314 0.242
8 2.31285 7.90299 19.5202 0.4440672 1.51737408 3.7478784 0.354667 0.575443 0.925672 0.192
9 3.25927 5.07819 15.1546 0.85392874 1.33048578 3.9705052 0.53194 0.703916 0.990093 0.262
10 3.28088 7.04839 19.8235 0.81693912 1.75504911 4.9360515 0.437181 0.76528 0.974517 0.249 ];
dataUsed = [0,1,1,1,0,0,0,0,0,0,1];
dataVariances = [0,0.0109,0.0374,0.0745,0.0000737,0.0001531,0.0000872,.0787,.0304,.0792,0.0000154];
maxIter = 20;
convtest = 0.00001;

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

% 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(dataVariances) ~= mpotential)
    fprintf(1,'wrong number of relative data-variances  \n', size(dataVariances));
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;
% 
