% 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 = 8;
parNamesCell = {'  D11', '  D21', '  D31', '  Z1', '  D12', '  D22', '  D32', '  Z2'};  %Dic is demand for node i by class c
parInteger = [0,0,0,0,0,0,0,0]';  %all are reals
parInit = [1,1,2,10,1,1,1,30]';  %10 and 30 are fixed
parFit = [1,1,1,0,1,1,1,0]';  %the demands are fitted, the think times are not
pControlInd = [1,2];  %in the data matrix, col 1 and 2 are controlled parameters, the
% class sizes N1 and N2.

% Data-related:
TMAX = 56; %trials (the number of rows of data
mpotential = 22; %data elements per trial, the number of columns of data
%dataNames = [ '$N1', '$N2', '$T11', '$T21', '$T31', '$n11', '$n21', '$n31', '$U11', '$U21', '$U31', '$Thru1', '$T12', '$T22', '$T32', '$n12', '$n22', '$n32', '$U12', '$U22', '$U32','$Thru2']';

% simulated by round robin, fitted as processor sharing
data=[1 1 2.22148 3.65693 4.06217 0.112011 0.184159 0.203757 0.100954 0.151131 0.199631 0.050007 5.5158 11.5528 1.20809 0.114699 0.239981 0.0250737 0.103888 0.206663 0.0206497 0.0207713
4 1 2.93748 5.76302 7.1327 0.454983 0.892682 1.10456 0.310521 0.465694 0.617883 0.154741 7.23384 17.3272 2.17791 0.12723 0.303727 0.0386066 0.08858 0.17359 0.0175827 0.017655
7 1 3.54872 8.32575 12.2206 0.727751 1.71008 2.51063 0.409807 0.615278 0.824176 0.205394 8.44805 23.9372 3.79649 0.127716 0.360801 0.0574152 0.07584 0.149743 0.015077 0.015121
10 1 3.96906 10.898 18.9463 0.903986 2.48662 4.32435 0.455175 0.684766 0.913172 0.228364 9.3167 29.9521 5.90646 0.124115 0.395845 0.0785977 0.0668317 0.132673 0.0133693 0.0133433
13 1 4.20852 12.7868 27.5677 0.999063 3.04646 6.57754 0.475188 0.713822 0.955603 0.237707 9.63992 33.6403 8.66703 0.119115 0.40818 0.106671 0.0620117 0.12008 0.0125203 0.0121907
16 1 4.27702 14.4533 37.2339 1.03587 3.50353 9.03732 0.484521 0.726602 0.972448 0.242187 9.79781 36.4908 11.223 0.11187 0.415743 0.127921 0.0571483 0.114057 0.0114287 0.011399
19 1 4.26913 15.5081 47.8148 1.04663 3.79967 11.6995 0.4883 0.73577 0.981517 0.245378 9.7069 38.7247 14.1645 0.104651 0.421049 0.153802 0.0530033 0.108907 0.0108287 0.010781
1 5 3.21441 8.63548 4.37563 0.123638 0.327516 0.167347 0.0770993 0.113123 0.153611 0.0381723 8.02701 27.1324 1.26989 0.604826 2.04607 0.0959363 0.37582 0.75185 0.0758613 0.0752277
4 5 3.91921 13.5587 6.67044 0.460322 1.58708 0.780976 0.234044 0.351759 0.468345 0.117319 9.81895 41.6972 1.98669 0.586177 2.49212 0.118346 0.298273 0.5948 0.059348 0.060139
7 5 4.39851 19.0036 9.63723 0.715793 3.09218 1.56742 0.325991 0.487011 0.650248 0.162706 10.9586 58.5803 2.94492 0.533541 2.8612 0.144091 0.24342 0.48904 0.048737 0.048738
10 5 4.72311 24.2742 13.653 0.896968 4.60702 2.59733 0.380369 0.568172 0.762121 0.189863 11.7136 73.1138 4.27961 0.491594 3.06544 0.179464 0.209915 0.419507 0.041858 0.0420737
13 5 4.84182 29.6277 18.3567 1.00337 6.12579 3.79902 0.413516 0.619833 0.831795 0.207123 12.0664 88.0007 5.70431 0.44728 3.24146 0.210825 0.185118 0.37033 0.036791 0.0368363
16 5 5.00087 34.5253 23.5094 1.09324 7.56287 5.14855 0.438024 0.657968 0.875924 0.219396 12.4501 100.038 7.48722 0.414004 3.33426 0.249626 0.168027 0.33373 0.0333453 0.033408
19 5 5.08288 38.3852 30.0266 1.15618 8.731 6.84139 0.455426 0.683112 0.910999 0.227379 12.6737 109.628 9.41471 0.393186 3.38526 0.290611 0.156825 0.308723 0.0306947 0.0310347
1 9 3.8088 18.3307 4.40452 0.103285 0.501123 0.120775 0.054142 0.08243 0.109187 0.0273557 9.75005 58.7008 1.24478 0.878124 5.30333 0.112374 0.449703 0.907143 0.0899827 0.090312
4 9 4.18077 25.2117 6.01716 0.36646 2.2221 0.528116 0.174021 0.26451 0.351272 0.088022 10.6003 80.4022 1.70523 0.776937 5.89904 0.125137 0.365225 0.73257 0.0732683 0.07358
7 9 4.56242 32.1855 7.84879 0.586752 4.12265 1.01057 0.257392 0.385676 0.514389 0.128269 11.6715 103.364 2.27937 0.717406 6.30186 0.139957 0.309577 0.613097 0.0615383 0.0614543
10 9 4.66633 39.1661 10.0443 0.729453 6.13377 1.57279 0.31224 0.469119 0.626097 0.156826 11.8103 125.077 2.9669 0.625187 6.63064 0.156819 0.26481 0.53032 0.0529683 0.0528623
13 9 4.79617 45.9942 12.7855 0.846678 8.12696 2.26167 0.353277 0.529585 0.707523 0.176471 12.2208 146.237 3.80406 0.571488 6.84721 0.178162 0.234612 0.468873 0.04672 0.0467207
16 9 4.91532 52.5148 15.6965 0.947636 10.1073 3.01726 0.386404 0.57759 0.768919 0.192351 12.3446 166.794 4.72587 0.518927 7.02781 0.198169 0.20968 0.42296 0.042087 0.0419707
19 9 4.99047 58.4992 19.4301 1.02126 11.9706 3.9675 0.409988 0.613034 0.817232 0.203929 12.6654 183.817 5.83493 0.492591 7.12264 0.226992 0.194797 0.385933 0.0385157 0.0386273
1 13 4.04801 30.2084 4.44301 0.0828344 0.619383 0.091375 0.0412427 0.061704 0.0820067 0.0205443 10.0932 97.9642 1.22376 0.940985 9.15795 0.114701 0.466028 0.936613 0.0938547 0.0932753
4 13 4.29026 37.3581 5.65294 0.299514 2.60676 0.395673 0.140075 0.208688 0.280209 0.069899 10.9005 120.252 1.57143 0.870104 9.59276 0.125465 0.400287 0.791607 0.0800797 0.080062
7 13 4.4925 45.1238 6.99708 0.471884 4.74211 0.735823 0.210231 0.314973 0.421471 0.104945 11.2984 145.606 1.95416 0.779832 10.0319 0.134528 0.34569 0.687483 0.0690463 0.068712
10 13 4.61461 52.8086 8.47397 0.608308 6.9597 1.11602 0.263853 0.39517 0.526436 0.131656 11.6027 169.738 2.39795 0.707151 10.324 0.144832 0.303335 0.604997 0.060477 0.0607493
13 13 4.75322 60.0042 10.3049 0.727121 9.16572 1.57874 0.30588 0.457791 0.613176 0.152754 12.0462 194.972 2.91067 0.652345 10.5623 0.158003 0.271147 0.543593 0.053854 0.05436
16 13 4.76746 67.6706 12.1764 0.804084 11.4503 2.05243 0.337361 0.507431 0.674573 0.169005 11.96 219.233 3.49913 0.587951 10.7605 0.173263 0.244422 0.49172 0.0494 0.0492163
19 13 4.92887 74.8061 14.3664 0.899571 13.6667 2.6127 0.365416 0.548943 0.726836 0.182151 12.37 240.393 4.18385 0.558587 10.8893 0.189119 0.225583 0.45179 0.045566 0.0453737
1 17 4.01286 42.1307 4.47734 0.0662344 0.693762 0.0737569 0.033128 0.04939 0.065912 0.016567 10.1177 136.938 1.20415 0.964185 13.0506 0.114643 0.476418 0.95197 0.095423 0.0956087
4 17 4.24084 50.0092 5.46099 0.242977 2.86833 0.312719 0.114901 0.172419 0.230692 0.0574137 10.7153 162.924 1.46331 0.888987 13.5059 0.12173 0.41526 0.827 0.0833953 0.0829623
7 17 4.37728 58.2006 6.41915 0.386836 5.15625 0.570027 0.177308 0.267056 0.354084 0.0887383 11.0313 189.153 1.72772 0.807563 13.8684 0.126356 0.364818 0.733087 0.0728893 0.073398
10 17 4.49653 65.5447 7.52227 0.512487 7.48846 0.859129 0.228612 0.341676 0.455809 0.114103 11.2365 215.555 2.07326 0.738022 14.1507 0.136033 0.327217 0.659177 0.0656457 0.065494
13 17 4.63989 73.5346 8.98133 0.619936 9.84507 1.20204 0.267517 0.401553 0.535863 0.133749 11.6896 239.97 2.46969 0.702389 14.3591 0.147149 0.29953 0.598713 0.0592137 0.0597997
16 17 4.75146 81.0286 10.3047 0.717224 12.2168 1.55638 0.302319 0.452032 0.603791 0.15075 11.9479 264.215 2.84884 0.660897 14.5244 0.15747 0.275262 0.548387 0.054908 0.0551597
19 17 4.75448 89.0956 11.8205 0.780265 14.6399 1.9404 0.327673 0.492293 0.656788 0.163877 11.887 290.687 3.27657 0.599673 14.7154 0.165999 0.25189 0.507513 0.0504383 0.0507207
1 21 4.04234 54.8281 4.35014 0.0547998 0.750126 0.05923 0.0272107 0.041475 0.0538747 0.0135873 10.1048 178.345 1.17405 0.9664 17.0536 0.112234 0.480188 0.959463 0.0958013 0.0953363
4 21 4.2653 62.4098 5.19332 0.208723 3.04921 0.253893 0.0979447 0.146545 0.194548 0.048824 10.6595 203.494 1.39236 0.91313 17.4154 0.119517 0.429727 0.852547 0.0862333 0.085459
7 21 4.36398 69.9387 6.08483 0.338912 5.41554 0.472376 0.15541 0.232169 0.31138 0.0774913 10.9644 229.79 1.63009 0.847853 17.7007 0.12577 0.386653 0.76995 0.077078 0.0774787
10 21 4.44892 78.5971 6.95971 0.446758 7.84784 0.700421 0.200651 0.299883 0.4018 0.100451 11.2365 258.626 1.89931 0.783538 18.0013 0.132397 0.349558 0.69768 0.0697003 0.0694523
13 21 4.51841 86.4506 7.96884 0.540043 10.3121 0.950802 0.238587 0.357354 0.476821 0.119689 11.3906 283.452 2.16758 0.732073 18.1972 0.139462 0.32071 0.64346 0.0643097 0.064463
16 21 4.56997 94.575 9.01222 0.621636 12.8048 1.22059 0.271005 0.40605 0.540693 0.135346 11.4159 309.796 2.4496 0.676326 18.403 0.145299 0.296027 0.59438 0.0594013 0.0591727
19 21 4.63487 102.65 10.3711 0.69099 15.2732 1.54662 0.297208 0.44681 0.595567 0.149093 11.657 336.844 2.81826 0.638007 18.5697 0.154858 0.274122 0.55191 0.0550757 0.0547733
1 25 3.98204 67.4591 4.51859 0.0461512 0.786941 0.0523966 0.0229307 0.035188 0.0469293 0.011407 10.0454 218.516 1.17154 0.96429 21.0212 0.112735 0.479648 0.966117 0.0964253 0.0963253
4 25 4.17406 74.4985 5.11236 0.178375 3.17218 0.219405 0.0855633 0.127462 0.172243 0.0428417 10.502 245.117 1.34524 0.916763 21.3428 0.117819 0.434368 0.873527 0.087229 0.0874703
7 25 4.31057 83.0497 5.79343 0.292455 5.63965 0.392333 0.136019 0.2043 0.270923 0.0679083 10.7371 271.649 1.52908 0.854879 21.6245 0.122087 0.398107 0.796557 0.0797857 0.0799583
10 25 4.38083 91.2353 6.55472 0.391081 8.13294 0.586872 0.178311 0.268808 0.357317 0.089145 10.9174 298.398 1.75161 0.804929 21.8504 0.12882 0.366992 0.732177 0.0734647 0.0736927
13 25 4.44905 99.3801 7.28997 0.476371 10.669 0.780262 0.214399 0.32256 0.426151 0.107245 11.0734 325.385 1.95884 0.751344 22.0815 0.132653 0.338827 0.678963 0.0675867 0.067766
16 25 4.55969 107.638 8.25952 0.559941 13.2025 1.01444 0.245311 0.369466 0.490151 0.12239 11.4223 351.617 2.23054 0.723651 22.2333 0.14068 0.315932 0.63092 0.063111 0.0632473
19 25 4.65982 115.758 9.27745 0.632009 15.7466 1.26008 0.272414 0.408946 0.543531 0.13585 11.5074 379.076 2.48572 0.679291 22.4075 0.146974 0.294307 0.5913 0.0589763 0.058952
1 29 3.99425 78.0179 4.39392 0.0411709 0.807778 0.0462563 0.020404 0.030953 0.0416013 0.0104413 10.2043 257.704 1.15249 0.98843 24.9911 0.111644 0.485133 0.969543 0.0965447 0.096799
4 29 4.09137 87.027 5.04622 0.154137 3.27391 0.191256 0.0752373 0.113433 0.152612 0.037881 10.1545 286.581 1.32186 0.89907 25.3303 0.117126 0.440818 0.88768 0.088754 0.0883517
7 29 4.28609 94.8456 5.66739 0.262169 5.77871 0.346751 0.122539 0.181839 0.245367 0.061213 10.7768 313.265 1.49112 0.879744 25.5655 0.121299 0.408348 0.816833 0.0812957 0.0811033
10 29 4.32801 103.028 6.32414 0.35093 8.33273 0.510126 0.162443 0.242355 0.323509 0.0806797 10.7995 342.142 1.66952 0.814517 25.7852 0.126247 0.375523 0.757393 0.075629 0.0757193
13 29 4.4124 111.768 6.97146 0.430006 10.9163 0.680117 0.195273 0.292874 0.390133 0.0973823 11.0148 364.368 1.85561 0.784401 25.9588 0.132201 0.354895 0.710287 0.071121 0.0709283
16 29 4.49356 120.212 7.65505 0.504734 13.5111 0.859278 0.223757 0.337866 0.447388 0.112527 11.3108 393.754 2.05821 0.751723 26.1306 0.136116 0.332692 0.66262 0.0662507 0.066105
19 29 4.64845 127.941 8.66994 0.584274 16.0631 1.09253 0.251996 0.37683 0.505204 0.125972 11.7274 419.774 2.31289 0.737302 26.2387 0.144749 0.315018 0.624857 0.062609 0.062747];


%the following variances were taken from the confidence intervals reported
%by the simulation (square of half conf int/1.96)
dataVariances = [0,0,0.1425448, 8.7567680, 0.0703873, 0.0010998, 0.0031497, 0.0005271, 0.0000260, 0.0000128, 0.0001148, 0.0000023, 0.2061901, 10.0062474, 0.0250156, 0.0019252, 0.0093971, 0.0002041, 0.0001041, 0.0000128, 0.0000042, 0.0000034];

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

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

%parexact is the values used for generating data
%parexact = [2,5,3,10,10,4,1,30]';
% 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;
% 
